Inverse Mask Design and Correction for Electronic Design

ABSTRACT

Various implementations of the invention provide for the generation of “smooth” mask contours by inverse mask transmission derivation and by subsequently “smoothing” the derived mask contours by proximity correction.

RELATED APPLICATIONS

This application claims priority under 35 U.S.C. §119(e) to U.S.Provisional Patent Application No. 61/154,271 entitled “Extreme OpticalProcess Correction,” filed on Feb. 20, 2010, and naming Yuri Granik etal. as inventors, and is a continuation in part of U.S. application Ser.No. 12/416,016 entitled “Calculation System for Inverse Masks,” filedMar. 31, 2009, and naming Yuri Granik as inventor, which in turn claimsthe benefit of U.S. Provisional Patent Application No. 61/041,197, filedMar. 31, 2008 and is itself a continuation in part of U.S. applicationSer. No. 12/359,174, filed Jan. 23, 2009, which in turn claims thebenefit of U.S. Provisional Patent Application No. 60/792,476 filed Apr.14, 2006, and is a continuation in part of U.S. patent application Ser.No. 11/364,802 filed Feb. 28, 2006, which in turn claims the benefit ofU.S. Provisional Patent Application No. 60/657,260 filed Feb. 28, 2005;U.S. Provisional Patent Application No. 60/658,278, filed Mar. 2, 2005;and U.S. Provisional Patent Application No. 60/722,840 filed Sep. 30,2005, which applications are all incorporated entirely herein byreference.

FIELD OF THE INVENTION

The invention relates to the field of integrated circuit design andmanufacturing. More particularly, various implementations of theinvention are applicable to jointly calibrating models useful tosimulate optical lithographic masks.

BACKGROUND OF THE INVENTION

Electronic circuits, such as integrated microcircuits, are used in avariety of products, from automobiles to microwaves to personalcomputers. Designing and fabricating microcircuit devices typicallyinvolves many steps, sometimes referred to as the “design flow.” Theparticular steps of a design flow often are dependent upon the type ofmicrocircuit, its complexity, the design team, and the microcircuitfabricator or foundry that will manufacture the microcircuit. Typically,software and hardware “tools” verify the design at various stages of thedesign flow by running software simulators and/or hardware emulators.These steps aid in the discovery of errors in the design, and allow thedesigners and engineers to correct or otherwise improve the design.These various microcircuits are often referred to as integrated circuits(“IC 's”).

Several steps are common to most design flows. Initially, a design maytypically start at a high level of abstraction, by a designer creating aspecification that describes particular desired functionality. Thisspecification, often implemented by a programming language, such as, forexample, the C or C++ programming language, describes at a high levelthe desired behavior of the device. Designers will then typically takethis specification for the design and create a logical design, oftenimplemented in a netlist, through a synthesis process. The logicaldesign describes the individual components of the design, and also mayhave different level of abstraction, such as, for example the gate levelor the register level.

A register transfer level (“RTL”) design, often implemented by ahardware description language (“HDL”) such as Verilog, SystemVerilog, orVery High speed hardware description language (“VHDL”), describes theoperation of the device by defining the flow of signals or the transferof data between various hardware components within the design. Moreparticularly, a register transfer level design describes theinterconnection and exchange of signals between hardware registers andthe logical operations that are performed on those signals.

Typically, a register transfer level design is first synthesized fromthe specification, followed by a gate level design being synthesizedfrom the register transfer level design. Gate level designs describe theactual physical components such as transistors, capacitors, andresistors as well as the interconnections between these physicalcomponents. Often, gate level designs are also implemented by a netlist,such as, for example, a mapped netlist. Lastly, the gate-level design istaken and another transformation is carried out. First by place androute tools that arrange the components described by the gate-levelnetlist and route connections between the arranged components; andsecond, by layout tools that generate a layout description having layout“shapes” that may then used to fabricate the electronic device, throughfor example, an optical lithographic process.

Integrated circuit layout descriptions can be provided in many differentformats. The Graphic Data System II (“GDSII”) format is popular fortransferring and archiving two-dimensional graphical IC layout data.Among other features, it contains a hierarchy of structures, eachstructure containing layout elements (e.g., polygons, paths orpoly-lines, circles and textboxes). Other formats include an open sourceformat named Open Access, Milkyway by Synopsys, Inc., EDDM by MentorGraphics, Inc., and the more recent Open Artwork System InterchangeStandard (OASIS) proposed by Semiconductor Equipment and MaterialsInternational (“SEMI”). These various industry formats are used todefine the geometrical information in integrated circuit layout designsthat are employed to manufacture integrated circuits. Once themicrocircuit device design is finalized, the layout portion of thedesign can be used by fabrication tools to manufacturer the device usinga photolithographic process.

There are many different fabrication processes for manufacturing acircuit, but most processes include a series of steps that depositlayers of different materials on a substrate, expose specific portionsof each layer to radiation, and then etch the exposed (or non-exposed)portions of the layer away. For example, a simple semiconductor devicecomponent could be manufactured by the following steps. First, apositive type epitaxial layer is grown on a silicon substrate throughchemical vapor deposition. Next, a nitride layer is deposited over theepitaxial layer. Then specific areas of the nitride layer are exposed toradiation, and the exposed areas are etched away, leaving behind exposedareas on the epitaxial layer, (i.e., areas no longer covered by thenitride layer). The exposed areas then are subjected to a diffusion orion implantation process, causing dopants, for example phosphorus, toenter the exposed epitaxial layer and form charged wells. This processof depositing layers of material on the substrate or subsequent materiallayers, and then exposing specific patterns to radiation, etching, anddopants or other diffusion materials, is repeated a number of times,allowing the different physical layers of the circuit to bemanufactured.

Each time that a layer of material is exposed to radiation, a mask mustbe created to expose only the desired areas to the radiation, and toprotect the other areas from exposure. The mask is created from circuitlayout data. That is, the geometric elements described in layout designdata define the relative locations or areas of the circuit device thatwill be exposed to radiation through the mask. A mask or reticle writingtool is used to create the mask based upon the layout design data, afterwhich the mask can be used in a photolithographic process. The imageembodied in the layout data is often referred to as the intended ortarget image or target contours, while the image created in the mask isgenerally referred to as the mask contours. Furthermore, the imagecreated on the substrate by employing the mask in a photolithographicprocess is often referred to as the printed image or printed contours.

As designers and manufacturers continue to increase the number ofcircuit components in a given area and/or shrink the size of circuitcomponents, the shapes reproduced on the substrate become smaller andare placed closer together. The feature sizes are often referred to bythe distance between features, conventionally called the “process step,”or the “node.” For example, one nod is the 32 nanometer (“nm”) node.This implies that adjacent features in the design, such as, for example,identical cells in a memory array, are 32 nanometers apart. As processsteps are continually scaled down, the corresponding reduction infeature size increases the difficulty of faithfully reproducing theimage intended by the layout design onto the substrate. Various commontechniques exist for mitigating these pattern dependant effects. Thesetechniques are commonly referred to as “resolution enhancementtechniques” or “RETs.” For example, optical process correction (“OPC,”)douple patterning, and “assist features” (i.e. features inserted into amask that are not intended to be replicated during the manufacturingprocess, but which increase the fidelity of the image) are a few commonresolution enhancement techniques (RET) that are often employed toprepare a physical layout designs for manufacturing.

Even though conventional resolution enhancement techniques, such as, forexample optical process correction, provide for an increase in imagefidelity, this increase in image fidelity may not be sufficient atfuture nodes. For example, conventional optical process correctionlimits a particular mask shape to having a one to one comparison to thecorrected mask shape. Additionally, the topology of a mask shapecorrected by conventional resolution enhancement techniques will have asimilar topology to the original mask shape. Furthermore, the inclusionor insertion of assist features into a corrected mask is time consumingand requires significant effort on the part of the designer tofacilitate.

SUMMARY OF THE INVENTION

Various implementations of the invention provide for the generation of“smooth” mask contours by inverse mask transmission derivation and bysubsequently “smoothing” the derived mask contours by proximitycorrection.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention will be described by way of illustrativeembodiments shown in the accompanying drawings in which like referencesdenote similar elements, and in which:

FIGS. 1A-1G illustrate conventional Nashold projections;

FIGS. 2A-2C illustrate local minima for the structure shown in FIG. 11A;

FIGS. 3A-3D illustrate local minima for a pattern of contact holes;

FIGS. 4A-4B illustrate a comparison between local variations andgradient descent optimizations in accordance with an embodiment of thepresent invention;

FIGS. 5A-5E illustrate gradient descent solutions after 1, 5, 10, 20,and 50 iterations in accordance with an embodiment of the presentinvention;

FIGS. 6A-6C illustrate the results of solving an inverse mask problemwith contour fidelity metric for a positive mask in accordance with anembodiment of the present invention;

FIGS. 7A-7C illustrate a method of local variations with contourfidelity and PSM masks in accordance with another embodiment of thepresent invention;

FIGS. 8A-8C illustrate contact holes inserted around main contacts;

FIGS. 9A-9C illustrate contact holes for strong PSM masks;

FIGS. 10A-10B illustrate a layout inversion for random logic and an SRAMcell;

FIGS. 11A-11D illustrate a deconvolution by a Wiener filter inaccordance with one embodiment of the present invention;

FIGS. 12A-12D illustrate a result of a least square unconstrainedoptimization in accordance with an embodiment of the present invention;

FIGS. 13A-13D illustrate a result of a constrained least squareoptimization in accordance with an embodiment of the present invention;

FIGS. 14A-14C illustrate a method of local variations in accordance withan embodiment of the present invention;

FIGS. 15A-15C illustrate a method of local variations with contourfidelity in accordance with the present invention;

FIG. 16 illustrates a result of a local variations algorithm for a PSMmask in accordance with an embodiment of the present invention;

FIGS. 17A-17C illustrate a solution of phase and assist features inaccordance with an embodiment of the present invention;

FIGS. 18A-18F illustrate an example of the present invention as appliedto contact holes;

FIG. 19 illustrates a representative computer system that can be used toimplement the present invention;

FIGS. 20A and 20B illustrate a series of steps used to calculate a masklayout pattern in accordance with one embodiment of the presentinvention;

FIG. 21 illustrates one embodiment of a test pattern used to select anideal intensity threshold in accordance with an embodiment of thepresent invention;

FIG. 22 illustrates one optimized mask data pattern created to producethe test pattern shown in FIG. 21;

FIG. 23 illustrates the image intensity measured across a cutline of asimulated image of the test pattern shown in FIG. 21 for a conventionalOPC corrected mask and an optimized mask pattern created in accordancewith an embodiment of the present invention;

FIGS. 24A and 24B illustrate two techniques for creating mask featuresto correspond to the optimized mask data;

FIGS. 25A-25E illustrate an objective function and line search strategy;

FIGS. 26A and 26B illustrate a sub-step evaluation process for a linesearch strategy;

FIGS. 27A and 27B illustrate layout patterns with and without adaptiveweight processing for the inverse mask calculation;

FIGS. 28A and 28B illustrate layout patterns with and without adaptiveweight processing for the inverse mask calculation;

FIGS. 29A-29C illustrate mask inversion processing for target contacthole arrays;

DETAILED DESCRIPTION OF ILLUSTRATIVE IMPLEMENTATIONS

The operations of the disclosed implementations may be described hereinin a particular sequential order. However, it should be understood thatthis manner of description encompasses rearrangements, unless aparticular ordering is required by specific language set forth below.For example, operations described sequentially may in some cases berearranged or performed concurrently. Moreover, for the sake ofsimplicity, the illustrated flow charts and block diagrams typically donot show the various ways in which particular methods can be used inconjunction with other methods.

It should also be noted that the detailed description sometimes usesterms like “determine” to describe the disclosed methods. Such terms areoften high-level abstractions of the actual operations that areperformed. The actual operations that correspond to these terms willoften vary depending on the particular implementation, and will bereadily discernible by one of ordinary skill in the art.

Furthermore, in various implementations of the invention, a mathematicalmodel may be employed to represent an electronic device. With someimplementations, a model describing the connectivity of the device, suchas for example a netlist, is employed. Those of skill in the art willappreciate that the models, even mathematical models represent realworld device designs and real world physical devices. Accordingly,manipulation of the model, even manipulation of the model when stored ona computer readable medium, results in a different device design. Moreparticularly, manipulation of the model results in a transformation ofthe corresponding physical design and any physical device rendered ormanufactured by the device design. Additionally, those of skill in theart can appreciate that during many electronic design and verificationprocesses, the response of a devices design to various signals or inputsis simulated. This simulated response corresponds to the actual physicalresponse the device being modeled would have to these various signals orinputs.

Some of the methods described herein can be implemented by softwarestored on a computer readable storage medium, or executed on a computer.Accordingly, some of the disclosed methods may be implemented as part ofa computer implemented electronic design automation (EDA) tool. Theselected methods could be executed on a single computer or a computernetworked with another computer or computers. For clarity, only thoseaspects of the software germane to these disclosed methods aredescribed; product details well known in the art are omitted.

INTRODUCTION

As will be explained in further detail below, the present invention is amethod and apparatus for calculating a mask pattern that will print adesired layout or portion thereof on a wafer. FIG. 19 illustrates arepresentative computer system that can be used to calculate a masklayout pattern in accordance with one embodiment of the presentinvention. A computer system 50, including one or more programmableprocessors, receives a set of executable instructions on acomputer-readable media 52 such as a CD, DVD, or from a communicationlink such as a wired or wireless communication network such as theInternet. The computer system 50 executes the instructions to read allor a portion of a desired layout file from a database 54 or otherstorage media. The computer system 50 then calculates data for a masklayout by dividing a mask layout into a number of discrete pixels. Thecomputer system determines the transmission characteristic of each ofthe pixels so that the result on a wafer will substantially match thepattern defined in the desired layout file. After calculating the maskpixel transmission characteristics, the mask pixel data is used by amask writer 56 in order to produce one or more corresponding masks orreticles 58. In another embodiment of the invention, the computer system50 transmits the desired layout file or a portion thereof over acommunication link 60 such as the Internet, etc., to one or moreremotely located computers 62 that perform the mask data calculations.

FIGS. 20A-20B illustrate one sequence of steps used to calculate themask pixel data in accordance with an embodiment of the presentinvention. Although the steps are discussed in a particular order, itwill be appreciated by those skilled in the art that the steps may beperformed in a different order while still obtaining the functionalitydescribed.

Beginning at 100, a computer system obtains all or a portion of a layoutfile that defines a desired pattern of features to be created on asemiconductor wafer. At 102, the computer system divides the desiredlayout into a number of frames. In one embodiment, the frames formadjacent, partially overlapping areas in the layout. Each frame, forexample, may occupy an area of 5×5 microns. The size of each frame maydepend of the amount of memory available and the processing speed of thecomputer system being used.

At 104, the computer system begins processing each of the frames. At106, the computer defines a blank area of mask data that is pixilated.At 108, the computer defines a corresponding set of optimal mask data.In one embodiment, the optimal mask data defines a corresponding set ofpixels whose transmission characteristics are defined by the desiredlayout data. For example, each optimal mask data pixel in an area thatcorresponds to a wafer feature may have a transmission characteristic of0 (e.g. opaque) and each optimal mask data pixel that corresponds to anarea of no wafer feature may have a transmission characteristic of 1(e.g. clear). In some embodiments, it may be desirable to change thedata for the optimal mask data from that defined strictly by the desiredlayout pattern. For example, the corners of features may be rounded orotherwise modified to reflect what is practical to print on a wafer. Inaddition or alternatively, pixel transmission characteristics may bechanged from a binary 0/1 value to a grayscale value, to positive andnegative values (representing phase shifters) or to complex values(representing partial phase shifters).

In some embodiments the optimal mask data may also include or bemodified by a weighting function. The weighting function allows a userto determine how close the solution for a given pixel should be to thetransmission characteristic defined by the corresponding pixel in theoptimal mask data. The weighting function may be a number selectedbetween 0 and 1 that is defined for each pixel in the optimal mask data.

At 110, an objective function is defined that relates a simulation ofthe image intensity on wafer to the pixel transmission characteristicsof the mask data and the optics of the photolithographic printingsystem. The objective function may be defined for each frame of maskdata or the same objective function may be used for more than one frameof mask data. Typically, the objective function is defined so that thevalue of the objective is minimized with the best possible mask, howeverother possibilities could be used.

In one embodiment of the present invention, one or more penaltyfunctions are also defined for the objective function. The penaltyfunctions operate to steer the optimization routine described below to asolution that can be or is desired to be manufactured on a mask. Forexample, it may be that the objective function has a number of possiblesolutions of which only one can actually be made as a physical mask.Therefore, the penalty functions operate to steer the solution in adirection selected by the programmer to achieve a mask that ismanufacturable. Penalty functions can be defined that promote variousstyles of resolution enhancement techniques such as: assist features,phase shifters, partial phase shifters, masks having features withgrayscale transmission values or multiple transmission values,attenuated phase shifters, combinations of these techniques or the like.

For example, a particular mask to be made may allow each pixel to haveone of three possible values with a transmission characteristic of 0(opaque)+1 (clear) or −1 (clear with phase shift). By including apenalty function in the objective function prior to optimization, thesolution is steered to a solution that can be manufactured as this typeof mask.

An example of a penalty function is α₄∥(m+e)m(m−e)∥₂ ², where e is a onevector, as set forth in Equation 57 described in the “Fast Pixel-BasedMask Optimization for Inverse Lithography” paper below. In oneembodiment, the penalty functions are defined as polynomials havingzeroes at desired pixel transmission characteristics. In anotherembodiment, the penalty functions can represent logical operations. Forexample, if the area of a wafer is too dark, the corresponding pixels inthe mask data can be made all bright or clear. This in combination withother mask constraints has the effect of adding subresolution assistfeatures to the mask data.

At 112, the objective function, including penalty functions, for theframe is optimized. In one embodiment the optimized solution is foundusing a gradient descent. If the objective function is selected to havethe form described by Equation 57, its gradient can be mathematicallycomputed using convolution or cross-correlation, which is efficient toimplement on a computer. The result of the optimization is a calculatedtransmission characteristic for each pixel in the mask data for theframe.

At 114, it is determined if each frame has been analyzed. If not,processing returns to 104 and the next frame is processed. If all framesare processed, the mask pixel data for each of the frames is combined at116 to define the pixel data for one or more masks. The mask data isthen ready to be delivered to a mask writer in order to manufacture thecorresponding masks.

More mathematical detail of the method for computing the mask pixeltransmission characteristics is described U.S. Patent Application No.60/722,840 filed Sep. 30, 2005 and incorporated by reference herein, aswell as in the paper “Fast Pixel-Based Mask Optimization for InverseLithography” by Yuri Granik of Mentor Graphics Corporation, 1001 RidderPark Drive, San Jose, Calif. 95131, reproduced below (with slightedits).

The direct problem of optical microlithography is to print features on awafer under given mask, imaging system, and process characteristics. Thegoal of inverse problems is to find the best mask and/or imaging systemand/or process to print the given wafer features. In this study weproposed strict formalization and fast solution methods of inverse maskproblems. We stated inverse mask problems (or “layout inversion”problems) as non-linear, constrained minimization problems over domainof mask pixels. We considered linear, quadratic, and non-linearformulations of the objective function. The linear problem is solved byan enhanced version of the Nashold projections. The quadratic problem isaddressed by eigenvalue decompositions and quadratic programmingmethods. The general non-linear formulation is solved by the localvariations and gradient descent methods. We showed that the gradient ofthe objective function can be calculated analytically throughconvolutions. This is the main practical result because it enableslayout inversion on large scale in order of M log M operations for Mpixels.

The layout inversion goal appears to be similar or even the same asfound in Optical Proximity Correction (OPC) or Resolution EnhancementTechniques (RET). However, we would like to establish the inverse maskproblem as a mathematical problem being narrowly formulated, thoroughlyformalized, and strictly solvable, thus differentiating it from theengineering techniques to correct (“C” in OPC) or to enhance (“E” inRET) the mask. Narrow formulation helps to focus on the fundamentalproperties of the problem. Thorough formalization gives opportunity tocompare and advance solution techniques. Discussion of solvabilityestablishes existence and uniqueness of solutions, and guidesformulation of stopping criteria and accuracy of the numericalalgorithms.

The results of pixel-based inversions can be realized by the opticalmaskless lithography (OML) [31]. It controls pixels of 30×30 nm (inwafer scale) with 64 gray levels. The mask pixels can be negative toachieve phase-shifting.

Strict formulations of the inverse problems, relevant to themicrolithography applications, first appear in pioneering studies of B.E. A. Saleh and his students S. Sayegh and K. Nashold. In [32], Sayeghdifferentiates image restoration from the image design (a.k.a. imagesynthesis). In both, the image is given and the object (mask) has to befound. However, in image restoration, it is guaranteed that the image isachieved by some object. In image design the image may not be achievableby any object, so that we have to target the image as close as possibleto the desired ideal image. The difference is analogical to solving fora function zero (image restoration) and minimizing a function (imagedesign). Sayegh states the image design problem as an optimizationproblem of minimizing the threshold fidelity error F_(c) in trying toachieve given threshold θ at the boundary C of the target image ([32],p. 86):

$\begin{matrix}{{{F_{C}\lbrack {m( {x,y} )} \rbrack} = {\oint\limits_{C}{( {{I( {x,y} )} - \theta} )^{n}{{l\min}}}}},} & (1)\end{matrix}$

where n=2 and n=4 options are explored; I(x, y) is image from the maskm(x, y); x, y are image and mask planar coordinates. Optical distortionsare modeled by the linear system of convolution with the point-spreadfunction h(x, y), so that

I(x,y)=h(x,y)*m(x,y),  (2)

and for the binary mask (m(x,y)=0 or m(x,y)=1). Sayegh proposesalgorithm of one at a time “pixel flipping”. Mask is discretized, andthen pixel values 0 and 1 are tried. If the error (1) decreases, thenthe pixel value is accepted, otherwise it is rejected, and we try thenext pixel.

Nashold [22] considers a bandlimiting operator in the place of thepoint-spread function (2). Such formulation facilitates application ofthe alternate projection techniques, widely used in image processing forthe reconstruction and is usually referenced as Gerchberg-Saxton phaseretrieval algorithm [7]. In Nashold formulation, one searches for acomplex valued mask that is bandlimited to the support of thepoint-spread function, and also delivers images that are above thethreshold in bright areas B and below the threshold in dark areas D ofthe target:

x,yεB:I(x,y)>θ,

x,yεD:I(x,y)<θ  (4)

Both studies [32] and [22] advanced solution of inverse problems for thelinear optics. However, the partially coherent optics ofmicrolithography is not a linear but a bilinear system [29], so thatinstead of (2) the following holds:

I(x,y)=∫∫∫∫q(x−x ₁ ,x−x ₂ ,y−y ₁ ,y−y ₂)m(x ₁ ,y ₁)m*(x ₂ ,y ₂)dx ₁ dx ₂dy ₁ dy ₂,  (5)

where q is a 4D kernel of the system. While the pixel flipping [32] isalso applicable to the bilinear systems, Nashold technique relies on thelinearity. To get around this limitation, Pati and Kailath [25] proposeto approximate bilinear operator by one coherent kernel h, thepossibility that follows from Gamo results [6]:

I(x,y)≈λ|h(x,y)*m(x,y)|²,  (6)

where constant λ is the largest eigenvalue of q, and his thecorrespondent eigenfunction. With this the system becomes linear in thecomplex amplitude A of the electrical field

A(x,y)=√{square root over (λ)}h(x,y)*m(x,y).  (7)

Because of this and because h is bandlimited, the Nashold technique isapplicable.

Liu and Zakhor [19, 18] advanced along the lines started by the directalgorithm [32]. In [19] they introduced optimization objective as aEuclidean distance ∥•∥₂ between the target I_(ideal) and the actualwafer images

F _(I) [m(x,y)]=∥I(x,y)−I _(ideal)(x,y)∥₂→min.  (8)

This was later used in (1) as image fidelity error in sourceoptimization. In addition to the image fidelity, the study [18]optimized image slopes in the vicinity of the target contour C :

$\begin{matrix}{{{F_{S}\lbrack {m( {x,y} )} \rbrack} = {{\oint\limits_{C - ɛ}{{I( {x,y} )}{l}}} - {\oint\limits_{C + ɛ}{{I( {x,y} )}{{l\min}}}}}},} & (9)\end{matrix}$

where C+ε is a sized-up and C−ε s is a sized down contour C ; ε is asmall bias. This objective has to be combined with the requirement forthe mask to be a passive optical element m(x,y)m*(x,y)≦1 or, usinginfinity norm ∥.∥_(∞)=max|.|, we can express this as

∥m(x,y)∥_(∞)≦1.  (10)

In case of the incoherent illumination

I(x,y)=h(x,y)²*(m(x,y)m*(x,y))  (12)

the discrete version of (9,10) is a linear programming (LP) problem forthe square amplitude p_(i)=m_(i)m_(i)* of the mask pixels, and wasaddressed by the “branch and bound” algorithm. When partially coherentoptics (4) is considered, the problem is complicated by the interactionsm_(i)m_(j)* between pixels and becomes a quadratic programming (QP)problem. Liu [18] applied simulated annealing to solve it. Consequently,Liu and Zakhor made important contributions to the understanding of theproblem. They showed that it belongs to the class of the constrainedoptimization problems and should be addressed as such. Reduction to LPis possible; however, the leanest relevant to microlithography andrigorous formulation must account for the partial coherence, so that theproblem is intrinsically not simpler than QP. New solution methods, moresophisticated than the “pixel flipping,” have also been introduced.

The first pixel-based pattern optimization software package wasdeveloped by Y.-H. Oh, J-C Lee, and S. Lim [24], and called OPERA, whichstands for “Optical Proximity Effect Reducing Algorithm”. Theoptimization objective is loosely defined as “the difference between theaerial image and the goal image,” so we assume that some variant of (7)is optimized. The solution method is a random “pixel flipping,” whichwas first tried in [32]. Despite the simplicity of this algorithm, itcan be made adequately efficient for small areas if image intensity canbe quickly calculated when one pixel is flipped. The drawback is thatpixel flipping can easily get stuck in the local minima, especially forPSM optimizations. In addition, the resulting patterns often havenumerous disjoined pixels, so they have to be smoothed, or otherwisepost-processed, to be manufacturable [23]. Despite these drawbacks, ithas been experimentally proven in [17] that the resulting masks can bemanufactured and indeed improve image resolution.

The study [28] of Rosenbluth, A., et al., considered mask optimizationas a part of the combined source/mask inverse problem. Rosenbluthindicates important fundamental properties of inverse mask problems,such as non-convexity, which causes multiple local minima. The solutionalgorithm is designed to avoid local minima and is presented as anelaborate plan of sequentially solving several intermediate problems.

Inspired by the Rosenbluth paper and based on his dissertation and theSOCS decomposition [2], Socha delineated the interference mappingtechnique [34] to optimize contact hole patterns. The objective is tomaximize sum of the electrical fields A in the centers (x_(k), y_(k)) ofthe contacts k=1 . . . N:

$\begin{matrix}{{F_{B}\lbrack {m( {x,y} )} \rbrack} = {- {\sum\limits_{k}{{{A( {x_{k},y_{k}} )}\min}.}}}} & (13)\end{matrix}$

Here we have to guess the correct sign for each A(x_(k), y_(k)), becausethe beneficial amplitude is either a large positive or a large negativenumber ([34] uses all positive numbers, so that the larger A thebetter). When kernel h of (7) is real (which is true for the unaberratedclear pupil), A and F_(B) are also real-valued under approximation (7)and for the real mask m. By substituting (7) into (13), we get

$\begin{matrix}\begin{matrix}{{- {\sum\limits_{k}{A( {x_{k},y_{k}} )}}} = {{- {\sum\limits_{k}( {h*m} )}}_{{x = x_{k}},{y = y_{k}}}}} \\{= {{\sum\limits_{k}{( {h*m} ) \cdot {\delta ( {{x - x_{k}},{x - x_{k}}} )}}} =}} \\{{= {{- ( {h*m} )} \cdot ( {\sum\limits_{k}{\delta ( {{x - x_{k}},{x - x_{k}}} )}} )}},}\end{matrix} & (14)\end{matrix}$

where the dot denotes an inner product ƒ·g=∫∫fgdxdy. Using the followingrelationship between the inner product, convolution*, andcross-correlation ∘ of real functions

(ƒ*g)·p=ƒ·(g∘p),  (15)

we can simplify (14) to

$\begin{matrix}{{{- {\sum\limits_{k}{A( {x_{k},y_{k}} )}}} = {{{- ( {h\mspace{11mu} \bullet {\sum\limits_{k}{\delta ( {{x - x_{k}},{x - x_{k}}} )}}} )} \cdot m} = {{- G_{b}} \cdot m}}},} & (16)\end{matrix}$

where function G_(I) is the interference map [34]. With (16) the problem(13) can be treated as LP with simple bounds (as defined in [8]) for themask pixel vector m={m_(i)}

−G _(b) ·m→min

−1≦m_(i)≦1.  (17)

In an innovative approach to the joined mask/source optimization byErdmann, A., et al. [4], the authors apply genetic algorithms (GA) tooptimize rectangular mask features and parametric source representation.GA can cope with complex non-linear objectives and multiple localminima.

Reduction to Linear and Thresholding Operators

The inverse mask problem can be reduced to a linear problem as it isshown above for IML or in [11]. This however requires substantialsimplifications. Perhaps richer and more interesting is modeling with alinear system and thresholding.

The linearization (7) can be augmented by the threshold operator tomodel the resist response. Inverse problems for such systems can besolved by Nashold projections [22]. Nashold projections belong to theclass of the image restoration techniques, rather than to the imageoptimizations, meaning that the method might not find the solution(because it does not exists at all), or in the case when it doesconverge, we cannot state that this solution is the best possible. Ithas been noted in [30] that the solutions strongly depend on initialguesses and do not deliver the best phase assignment unless thealgorithm is steered to it by a good initial guess. Moreover, if theinitial guess has all phases set to 0, then so has the solution.

Nashold projections are based on Gerchberg-Saxton [7] phase retrievalalgorithm. It updates the current mask iterate m^(k) via

m ^(k+1)=(P _(m) P _(s))m ^(k),  (31)

where P_(s) is a projection operator into the frequency support of thekernel h, and P_(m) is a projection operator that forces thethresholding (4). Gerchberg-Saxton iterations tend to stagnate. Fienap[5] proposed basic input-output (BIO) and hybrid input-output (HIO)variations that are less likely to be stuck in the local minima. Thesevariations can be generalized in the expression

m ^(k+1)=(P _(m) P _(s)+α(γ(P _(m) P _(s) −P _(s))−P _(m+I))m^(k),  (32)

where I is an identity operator; α=1,γ=0 for BIO, α=1,γ=1 for HIO, andα=0, γ=0 for the Gerchberg-Saxton algorithm.

We implemented operator P_(m) as a projection onto the ideal image

$\begin{matrix}{{{P_{m}m^{k}} = {\frac{m^{k}}{m^{k}}\sqrt{I_{ideal}}}},} & (33)\end{matrix}$

and P_(s) as a projection to the domain of the kernel h, i.e. P_(s)zeros out all frequencies of m which are high than the frequencies ofthe kernel h. The iterates (32) are very sensitive to the values of itsparameters and the shape of ideal image. We have found meaningfulsolutions only when the ideal image is smoothed. Otherwise the phasescome out “entangled,” i.e. the phase alternates along the lines as inFIG. 1E, right, instead of alternating between lines. We used Gaussiankernel with the diffusion length of 28 nm, which is slightly larger thanthe pixel size 20 nm in our examples. The behavior of iterates (32) isnot yet sufficiently understood [36], which complicates choice of α,γ.In our examples the convergence is achieved for α=0.9, γ=1 after T=5000iterations. When α=0, γ=0, which corresponds to the original Nasholdprojections (31), the iterations quickly stagnate converging to anon-printable mask. The runtime is proportional to T*M*log M, M- numberof pixels. The convergence is slow because T is large, so thatapplication to the large layout areas is problematic.

As shown in FIG. 1B, generalized Nashold projections (32) assignalternating phases to the main features and insert assists betweenlines. The lines widths are on-target, but line-ends are not corrected.The solution has good contrast. When projections stagnate, the phasesalternate along the lines. This “phase entanglement” is observedsometimes in the non-linear problems (considered in a section below)when their iterations start from the random pixel assignments.

Quadratic Problems

In the quadratic formulations of the inverse problems, the coherentlinearization (6) is not necessarily. We can directly use bilinearintegral (5). Our goal here is to construct objective function as is aquadratic form of mask pixels. We start with (8) and replace Euclideannorm (norm 2) with Manhattan norm (norm 1):

F _(I) [m(x,y)]=∥I(x,y)−I _(ideal)(x,y)∥₁→min.  (34)

The next step is to assume that the ideal image is sharp, 0 in darkregions and 1 in bright regions, so that I(x,y)≧I_(ideal)(x,y) in thedark regions and I(x,y)≦I_(ideal)(x,y) in the bright regions. This letsus to remove the module operation from the integral (34):

∥I(x,y)−I _(ideal)(x,y)∥₁ =∫∫|I−I _(ideal) |dxdy=∫∫w(x,y)(I(x,y)−I_(ideal)(x,y))dxdy,  (35)

where w(x, y) is 1 in dark regions and −1 in bright regions. Finally wecan ignore the constant term in (35), which leads to the objective

F _(w) [m(x,y)]=∫∫wI(x,y)→min.  (36)

The weighting function w can be generalized to have any positive valuein dark regions, any negative value in bright regions, and 0 in theregions which we choose to ignore. Proper choice of this function coversthe image slope objective (9), but not the threshold objective (1).Informally speaking, we seek to make bright regions as bright aspossible, and dark regions as dark as possible. Substituting (5) into(36), we get

∫∫wI(x,y)dxdy=∫∫∫∫Q(x ₁ y ₁ ,x ₂ ,y ₂)m*(x ₂ ,y ₂)dx ₁ dx ₂ dy ₁ dy₂,  (37)

Where

Q(x ₁ y ₁ ,x ₂ ,y ₂)=∫∫w(x,y)q(x−x ₁ ,x−x ₂ ,y−y ₁ ,y−y ₂)dxdy.  (38)

Discretization of (37) results to the following constrained QP

F _(w) [m]=m*Qm→min

∥m∥_(∞)≦1.  (39)

The complexity of this problem depends on the eigenvalues of matrix Q .When all eigenvalues are non-negative, then it is convex QP and anylocal minimizer is global. This is a very advantageous property, becausewe can use any of the numerous QP algorithms to find the global solutionand do not have to worry about local minima Moreover, it is well knownthat the convex QP can be solved in polynomial time. The next case iswhen all eigenvalues are non-positive, a concave QP. If we removeconstraints, the problem becomes unbounded, with no minimum and nosolutions. This means that the constraints play a decisive role: allsolutions, either local or global, end up at some vertex of the box∥m∥_(∞)≦1. In the worst case scenario, the solver has to visit allvertices to find the global solution, which means that the problem isNP-complete, i.e. it may take an exponential amount of time to arrive atthe global minima. The last case is an indefinite QP when both positiveand negative eigenvalues are present. This is the most complex andintractable case: an indefinite QP can have multiple minima, all lie onthe boundary.

We conjecture that the problem (39) belongs to the class of indefiniteQP. Consider case of the ideal coherent imaging, when Q is a diagonalmatrix. Vector w lies along its diagonal. This means that eigenvaluesμ₁, μ₂ . . . of Q are the same as components of the vector w, which arepositive for dark pixels and negative for bright pixels. If there is atleast one dark and one bright pixel, the problem is indefinite. Anotherconsideration is that if we assume that (39) is convex, then thestationary internal point m=0 where the gradient is zero

$\begin{matrix}{\frac{{F_{w}\lbrack m\rbrack}}{m} = {{2{Qm}} = 0}} & (40)\end{matrix}$

is the only solution, which is a trivial case of mask being dark. Thismeans that (39) is either has only a trivial solution, or it isnon-convex.

Related to (39) QP was considered by Rosenbluth [28]:

m*Q _(d) m→min

m*Q _(b) m≧b,  (41)

where Q_(d) and Q_(b) deliver average intensities in bright and darkregions correspondingly. The objective is to keep dark regions as darkas possible while maintaining average intensity not worse than somevalue b in bright areas. Using Lagrange multipliers, we can convert (41)to

m*(Q _(d) −λQ _(b))→min

∥m∥_(∞)≦1

λ≧0,  (42)

which is similar to (39).

Another metric of the complexity of (39) is number of the variables,i.e. the pixels in the area of interest. According to Gould [10], theproblems with order of 100 variables are small, more than 10³ are large,and more than 10⁵ are huge. Considering that the maskless lithographycan control transmission of the 30 nm by 30 nm pixel [31], the QP (39)is large for the areas larger than 1 um by 1 um, and is huge for theareas lager than 10 um by 10 um. This has important implications for thetype of the applicable numerical methods: in large problems we can usefactorizations of matrix Q, in huge problems factorizations areunrealistic.

For the large problems, when factorization is still feasible, a dramaticsimplification is possible by replacing the infinity norm by theEuclidean norm in the constraint of (39), which results in

F _(w) [m]=m*Qm→min

∥m∥₂≦1.  (43)

Here we search for the minimum inside a hyper-sphere versus a hyper-cubein (39). This seemingly minor fix carries the problem out of the classof NP-complete to P (the class of problems that can be solved inpolynomial time). It has been shown in [35] that we can find globalminima of (43) using linear algebra. This result served as a base forthe computational algorithm of “trust region” [13] which specificallyaddresses indefinite QP.

The problem (43) has the following physical meaning: we optimize thebalance of making bright regions as bright as possible and dark regionsas dark as possible while limiting light energy ∥m∥₂ ² coming throughthe mask. To solve this problem, we use procedures outlined in [35, 13].First we form Lagrangian function of (43)

L(m,λ)=m*Qm+λ(∥m∥ ²−1).  (44)

From here we deduce the first order necessary optimality conditions ofKarush-Kuhn-Tucker (or KKT conditions, [12]):

2(Q+λI)m=0

λ(∥m∥−1)=0

λ≧0

∥m∥≦1.  (45)

Using Sorensen [35], we can state what that (43) has a global solutionif and only if we can find such λ and m that (45) is satisfied and thematrix Q+λI is positive semidefinite or positively defined. Let us findthis solution.

First we notice that we have to choose λ large enough to compensate thesmallest (negative) eigenvalue of Q, i.e.

λ≧|μ₁|>0.  (46)

From the second condition in (45) we conclude that ∥m∥=1, that is thesolution lies on the surface of hyper-sphere and not inside it. The lastequation to be satisfied is the first one from (45). It has anon-trivial ∥m∥>0 solution only when the lagrange multiplier λ equals toa negative of one of the eigenvalues λ=−μ_(i). This condition and (46)has a unique solution λ=−μ₁, because other eigenvalues μ₂, μ₃, . . . areeither positive so that λ≧0 does not hold, or they are negative, butwith absolute value that is smaller than μ₁, so that λ≧|μ₁| not hold.

After we determined that λ=−μ₁, we can find m from 2(Q−μ₁I)m=0 as thecorresponding eigenvector m=v₁. This automatically satisfies ∥m∥=1,because all eigenvectors are normalized to have a unit length. Weconclude that (43) has a global solution which corresponds to thesmallest negative eigenvalue of Q .

As we have shown, the minimum eigenvalue of Q and its eigenvector playspecial role in the problem by defining the global minimum. However,other negative eigenvectors are also important, because it is easy tosee that any pair

λ=−μ_(i)>0

m=v_(i)  (47)

is a KKT point and as such defines a local minimum. The problem has asmany local minima as negative eigenvalues. We may also consider startingour numerical minimization from one of these “good” minima, because itis possible that a local minimum leads to a better solution in thehyper-cube than a global minimum of the spherical problem.

FIGS. 2A, 2B, and 2C shows three strongest local minima of the problem(39) for the structure of FIG. 1B. These local minima are pointing tothe beneficial interactions between layout features, suggestingalternating phase assignments. For example, the second solution suggeststhat the L-shape transmission should be chosen positive, while the combhas negative transmission, the dense vertical line of the comb haspositive transmission, and the second horizontal line has negativetransmission.

Results of the similar analysis for the case of the contact holes aredisplayed in FIG. 3A-3D. These results are stronger, and can be useddirectly in applications. The method “proposes” beneficial phases forthe contacts and position and phases of the assists. The mostinteresting solution is shown in the low right inset, where all contactshave well-defined transmissions, with 3 contacts positive and 4 contactsnegative. The advantages of this method comparing to IML [34] is thatthis method automatically finds the best phases of the contacts and isnot based of the coherent approximation.

FIG. 3B, 3C, 3D illustrates the first three local minima for QP on ahyper-sphere for the contact holes and process conditions from Socha[34]. The third solution has the clearest phase assignments and theposition of assists.

For the positive masks, in particular for the binary masks, theconstraint has to be tightened to ∥m−0.5∥_(∞)≦0.5. Then thecorrespondent to (39) problem is

F _(w) [m]=m*Qm→min

∥m−0.5∥_(∞)≦0.5.  (48)

This is also an indefinite QP and is NP-complete. Replacing hereinfinity norm with Euclidean norm, we get a simpler problem

m*Qm→min

∥Δm∥₂≦0.5

Δm=m−m ₀ ,m ₀={0.5,0.5, . . . , 0.5}.  (49)

The Lagrangian can be written as

L(m,λ)=m*Qm+λ(∥m−m ₀∥²−0.25).  (50)

The KKT point must be found from the following conditions

(Q+λI)Δm=−Qm ₀

λ(∥Δm∥ ²−0.25)=0

λ≧0

∥Δm∥≦0.5.  (51)

This is more complex problem than (45) because the first equation is nothomogeneous and the pairs λ=μ_(i), Δm =v_(i) are clearly not thesolutions. We can still apply the condition of the global minimumλ≧−μ₁>0 (Sorensen [35]). From the second condition we conclude that∥Δm∥²=0.25, meaning that all solutions lie on the hyper-sphere with thecenter at m₀. The case λ=−μ₁ is eliminated because the first equation isnot homogeneous, so that we have to consider only λ>−μ₁. Then Q+λI isnon-singular, we can invert it, and find the solution

Δm=−(Q+λI)⁻¹ Qm ₀.  (52)

The last step is to find the Lagrange multiplier λ that satisfy theconstraint ∥Δm∥²=0.25, that is we have to solve

∥(Q+λI)⁻¹ Qm ₀∥=0.5.  (53)

The norm on the right monotonically increases from 0 to infinity in theinterval −∞<λ<−μ₁, thus (53) has to have exactly one solution in thisinterval. The pair λ,Δm that solves (52-53) is a global solution of(49). We conjecture that there are fewer KKT points of local minima of(49) than in (45) (may be there are none), but this remains to be provenby analyzing behavior of the norm (53) when Lagrange multiplier isbetween negative eigenvalues. The solutions of (49) are supposed to showhow to insert assist features when all contacts have the same phases.

General Non-Linear Problems

Consider objective (8) of image fidelity error

F _(I) [m(x,y)]=∥I(x,y)−I _(ideal)(x,y)∥→min.  (54)

We can state this in different norms, Manhattan, infinity, Euclidean,etc. The simplest case is a Euclidean norm, because (54) becomes apolynomial of the forth degree (quartic polynomial) of mask pixels. Theobjective function is very smooth in this case, which ease applicationof the gradient-descent methods.

We can generalize (54) by introducing weighting w=w(x,y) to emphasizeimportant layout targets and consider smoothing in Sobolev norms as in[20]:

F _(w) [m(x,y)]²=∥√{square root over (w)}·(I−I _(ideal))∥₂ ²+α₁ ∥L ₁ m∥₂ ²+α₃ ∥m−m ₀∥₂ ²→min,  (55)

where L₁, L₂ are the operators of first and second derivatives,m₀=m₀(x,y) is some preferred mask configuration that we want to be closeto (for example, the target), and α₁, α₂, α₃ are smoothing weights. Thesolutions of (55) increase image fidelity; however, the numericalexperiments show that the contour fidelity of the images is notadequate. To address, we explicitly add (1) into (55):

$\begin{matrix}{{F_{wc}\lbrack m\rbrack}^{2} = {{{{\sqrt{w} \cdot ( {I - I_{ideal}} )}}_{2}^{2} + {\oint\limits_{C}{( {I - \theta} )^{n}{l}}} + {\alpha_{1}{{L_{1}m}}_{2}^{2}} + {\alpha_{2}{{L_{2}m}}_{2}^{2}} + {\alpha_{3}{{{m - m_{0}}}_{2}^{2}\min}{m}_{\infty}}} \leq 1}} & (56)\end{matrix}$

If the desired output is a two-, tri-, any other multi-level tone mask,we can add penalty for the masks with wrong transmissions. The simplestform of the penalty is a polynomial expression, so for example for thetri-tone Levenson-type masks with transmissions −1, 0, and 1, weconstruct the objective as

$\begin{matrix}{{{F_{wce}\lbrack m\rbrack}^{2} = {{{{\sqrt{w} \cdot ( {I - I_{ideal}} )}}_{2}^{2} + {\oint\limits_{C}{( {I - \theta} )^{n}{{l++}}\alpha_{1}{{L_{1}m}}_{2}^{2}}} + {\alpha_{2}{{L_{2}m}}_{2}^{2}} + {\alpha_{3}{{{m - m_{0}}}_{2}^{2}++}\alpha_{4}{{( {m + e} ){m( {m - e} )}}}_{2}^{2}{m}_{\infty}}} \leq {1\min}}},} & (57)\end{matrix}$

where e is a one-vector. Despite all the complications, the objectivefunction is still a polynomial of the mask pixels. To optimize for thefocus depth, the optimization of (57) can be conducted off-focus, as wassuggested in [16, 20]. After discretization, (55) becomes a non-linearprogramming problem with simple bounds.

We expect that this problem inherits property of having multiple minimafrom the corresponding simpler QP, though smoothing operators of (57)have to increase convexity of the objective. In the presence of multiplelocal minima the solution method and staring point are highlyconsequential: some solvers tend to converge to the “bad” localsolutions with disjoined masks pixels and entangled phases, othersbetter navigate solution space and chose smoother local minima. TheNewton-type algorithms, which rely on the information about secondderivatives, should be used with a caution, because in the presence ofconcavity in (57), the Newtonian direction may not be a descentdirection. The branch-and-bound global search techniques [18] are notthe right choice because they are not well-suited for the largemulti-dimensional optimization problems. It is also tempting to performnon-linear transformation of the variables to get rid of the constraintsand convert problem to the unconstrained case, for example by usingtransformation x_(i)=tanh(m_(i)) or m_(i)=sin(x_(i)) as in [26].

Solution Methods

The reasonable choices to solve (57) are descent algorithms withstarting points found from the analytical solutions of the related QP.We apply algorithms of local variations (“one variable at a time”),which is similar in spirit to the pixel flipping [32, 17], and also usea variation of the steepest descent by Frank and Wolfe [21] to solveconstrained optimization problems.

In the method of local variation, we chose the step Δ₁to compare threeexploratory transmissions for the pixel i: m_(i) ¹, m_(i) ¹+Δ₁, andm_(i) ¹−Δ₁. If one of these values violates constraints, then it ispulled back to the boundary. The best of these three values is accepted.We try all pixels, optionally in random exhaustive or circular order,until no further improvement is possible. Then we reduce step Δ₂<Δ₁ andrepeat the process until the step is deemed sufficiently small. Thisalgorithm is simple to implement. It naturally takes care of the simple(box) constraints and avoids the general problem of other moresophisticated techniques, which may converge prematurely to anon-stationary point. This algorithm calculates the objective functionnumerous times; however, the runtime cost of its exploratory calls isrelatively low with the electrical field caching (see the next section).Other algorithms may require fewer but more costly non-exploratorycalls. This makes method of local variation a legitimate tool in solvingthe problem, though descent methods that use convolution for thegradient calculations are faster.

Frank and Wolfe method is an iterative gradient descent algorithm tosolve constrain problems. At each step k we calculate the gradient∇F^(k) of the objective and then replace the non-linear objective withits linear approximation. This reduces the problem to LP with simplebounds:

∇F ^(k) ·m→min

∥m∥_(∞)1.  (59)

The solution of this m=l^(k) is used to determine the descent direction

p ^(k) =l ^(k) −m ^(k−1).  (60)

Then the line search is performed in the direction of p^(k) to minimizethe objective as a function one variable γε[0,1]:

F[m ^(k−1) +γp ^(k)]→min.  (61)

The mask M^(k)=m^(k−1)+γp^(k) is accepted as the next iterate. Theiterations continue until convergence criteria are met. Electrical fieldcaching helps to speedup line search and the gradient calculations ifnumerical differentiation is used.

FIGS. 4A and 4B show a comparison between (FIG. 4A) local variations andgradient descent optimizations (FIG. 4B) for the target as an isolatedcontact hole of 100 nm size, printed with quadrupole illumination. Themask pixels of 10 nm size are thresholded at 0.5 level (the dark areashave transmissions >0.5) and cover 2.56 by 2.65 um layout area. Thesimulation time of gradient descent with analytically calculatedgradient is about 2 seconds on a typical workstation. Solution withlocal variations is less regular than with the gradient descent, becausepixels are iterated randomly. Up to 8 assist rinds are discernable.

FIGS. 5A-5E show gradient descent mask iterations after 1, 5, 10, 20,and 50 iterations. The assist features become more and more complicatedas the descent iterations improve objective function. This is simulatedunder the same process conditions and a target layout as in FIG. 4B.

Gradient of Objective Function

The gradient descent algorithms require recalculation of the objectiveand its gradient at each iteration. The gradient of the objectivefunction can be calculated numerically or analytically. When theobjective is expressed in norm 2 as in (55), the derivatives can becalculated analytically, yielding efficient representation throughconvolutions.

Consider objective in the form of the weighted inner product(f,g)=∫∫wfgdxdy:

F _(w) ² [m]=∥√{square root over (w)}·(I−I _(ideal))∥²=(I−I _(ideal),I−I _(ideal)).  (63)

Small variations δm of the mask m cause the following changes in theobjective:

$\begin{matrix}\begin{matrix}{{F_{w}^{2}\lbrack {m + {\delta \; m}} \rbrack} = {( {{{I( {m + {\delta \; m}} )} - I_{ideal}},{{I( {m + {\delta \; m}} )} - I_{ideal}}} ) =}} \\{= {( {{{I(m)} + {\delta \; I} - I_{ideal}},{{I(m)} + {\delta \; I} - I_{ideal}}} ) \approx}} \\{\approx {{F_{w}^{2}\lbrack m\rbrack} + {2{( {{I - I_{ideal}},{\delta \; I}} ).}}}}\end{matrix} & (64)\end{matrix}$

Let us find δI =I(m+δm)−I(m). Using SOCS formulation (60), andneglecting O(δm²) terms, we get

$\begin{matrix}\begin{matrix}{{\delta \; I} = {{I( {m + {\delta \; m}} )} - {I(m)}}} \\{= {{\sum\limits_{i = 1}^{N}{{\lambda_{i}( {h_{i}*( {m + {\delta \; m}} )} )}( {h_{i}*( {m + {\delta \; m}} )} )^{*}}} -}} \\{{{\sum\limits_{i = 1}^{N}{\lambda_{i}( {h_{i}*m} )( {h_{i}*m} )^{*}}} \approx}} \\{\approx {\sum\limits_{i = 1}^{N}{\lambda_{i}\lbrack {{A_{i}( {h_{i}*\delta \; m} )}^{*} + {A_{i}^{*}( {h_{i}*\delta \; m} )}} \rbrack}}}\end{matrix} & (65)\end{matrix}$

where A_(i) is defined in (60). To use this in (64), we have to findscalar product of δI with ΔI=I−I_(ideal):

$\begin{matrix}\begin{matrix}{( {{\Delta \; I},{\delta \; I}} ) = {{\sum\limits_{i = 1}^{N}{\lambda_{i}\lbrack {( {{\Delta \; I},{A_{i}( {h_{i}*\delta \; m} )}^{*}} ) + ( {{\Delta \; I},{A_{i}^{*}( {h_{i}*\delta \; m} )}} )} \rbrack}} =}} \\{= {{\sum\limits_{i = 1}^{N}{\lambda_{i}\lbrack {( {{A_{i}\Delta \; I},{h_{i}^{*}*\delta \; m^{*}}} ) + ( {{A_{i}^{*}\Delta \; I},{h_{i}*\delta \; m}} )} \rbrack}} =}} \\{= {2{\sum\limits_{i = 1}^{N}{\lambda_{i}{{Re}( {{A_{i}^{*}\Delta \; I},{h_{i}*\delta \; m}} )}}}}}\end{matrix} & (66)\end{matrix}$

Using the following property of the weighted inner product

(f*g,h)=f·(g*∘wh)  (67)

we can convert (66) to the form

$\begin{matrix}\begin{matrix}{( {{\Delta \; I},{\delta \; I}} ) = {2{\sum\limits_{i = 1}^{N}{\lambda_{i}{{Re}( {{\delta \; m*h_{i}},{A_{i}^{*}\Delta \; I}} )}}}}} \\{= {2{\sum\limits_{i = 1}^{N}{\lambda_{i}{{Re}( {\delta \; {m \cdot ( {h_{i}^{*}\bullet \; w\; A_{i}^{*}\Delta \; I} )}} )}}}}}\end{matrix} & (68)\end{matrix}$

Substituting this into (64) gives us an analytical expression for thegradient of the objective

$\begin{matrix}{{{{F_{w}^{2}\lbrack {m + {\delta \; m}} \rbrack} - {F_{w}^{2}\lbrack m\rbrack}} \approx {{{\nabla F_{w}^{2}} \cdot \delta}\; m}}{{\nabla F_{w}^{2}} = {4{\sum\limits_{i = 1}^{N}{\lambda_{i}{{Re}( {h_{i}^{*}\bullet \; {wA}_{i}^{*}\Delta \; I} )}}}}}} & (69)\end{matrix}$

This formula let us calculate gradient of the objective throughcross-correlation or convolution as O(NM log(M)) FFT operation, which issignificantly faster than numerical differentiation with O(NM²)runtime.

Electrical Field Caching

The speed of the local variation algorithm critically depends on theability to quickly re-calculate image intensity when one or a few pixelschange. We use electrical field caching procedure to speed up thisprocess.

According to SOCS approximation [3], the image intensity is thefollowing sum of convolutions of kernels h_(i) (x, y) with the mask m(x,y):

$\begin{matrix}{{{I( {x,y} )} = {\sum\limits_{i = 1}^{N}{\lambda_{i}{A_{i}( {x,y} )}{A_{i}^{*}( {x,y} )}}}},{A_{i} = {{h_{i}( {x,y} )}*{{m( {x,y} )}.}}}} & ( {60A} )\end{matrix}$

Suppose that we know the electrical fields A_(i) ⁰ for the mask m⁰ andwant to calculate intensity for the slightly different mask m′. Then

A _(i) ′=A _(i) ⁰ +h _(i)*(m′−m ⁰)  (61A)

These convolutions can be quickly calculated by the directmultiplication, which is O(d·M·N) operation, where d is the number ofdifferent pixels between m⁰ and m′, M is pixel count of the kernels, andN is number of kernels. This may be faster than convolution by FFT.Constantly updating the cache A_(i) ⁰, we can quickly re-calculateintensities for small evolutionary mask changes.

The additivity of the electrical fields can also be exploited to speedupintensity calculations in the line search (61A). If the mask m^(k−1)delivers electrical fields A_(i) ^(k−1), and the mask p^(k) deliversB_(i) ^(k), then the intensity from the mask m=m^(k−1)+γp^(k) can bequickly calculated through its electrical fields A_(i):

A _(i) =A _(i) ^(k−1) +γB _(i) ^(k).  (62A)

This avoids convolutions of (60A) and reduces intensity calculation tomultiplication of the partial electrical fields A_(i).

Simulation Results

In FIGS. 6A-6C we show results of solving (56) with the contour fidelitymetric for the positive mask 0≦m≦1. The assist features can be seenaround main structures.

FIGS. 6A-6C show local variations method for the objective (57) withcontour fidelity. The contours are on target almost everywhere,including line ends. The image contrast is improved. Mask has rows ofassist features and serifs.

Next example demonstrates solutions when main features have the samephase and assist features can have phase shift, FIGS. 7A-7C. We observenegative transmission of the assists on the mask. The contrast along thecutline is improved in comparison to the ideal case (mask equal target).Contour fidelity is very good, the third inset. The last example iscontact holes, FIGS. 10A-10B and 11A-11D. The method is capable ofinserting assist contacts and deliver complex interferometric assistfeatures in PSM case.

FIGS. 7A-7C illustrate a method of local variations for the objective(57) with contour fidelity and PSM mask with assist features.

FIGS. 8A-8C illustrate contact holes example for the binary mask. Smallassist contact holes are inserted around main contacts. The imagecontrast is compared to the case when mask is the same as target. Thecontrast is improved significantly. Image contours are on target, thethird column.

FIGS. 9A-9C illustrate contact holes example for strong PSM mask.Resulting PSM mask has complex structure of assists holes, which arehard to separate from the main features. The contrast is even betterthan for the binary mask. Despite very complex mask, the contours are ontarget (lower right inset) and sidelobs do not print.

FIGS. 10A-10B demonstrate application of gradient descent to large pieceof the layout with contact holes. The target holes are shown in green.The resulting mask is thresholded at the transmission level 0.5. Themanufacturability of such complex masks is discussed in [33].

FIGS. 10A-10B illustrate examples of layout inversion for random logicand an SRAM cell.

We classified methods for solving inverse mask problems as linear,quadratic, and non-linear. We showed how to solve a quadratic problemfor the case of spherical constraint. Such analytical solutions can beused as a first step for solving non-linear problems. In the case of thecontacts, these solutions can be immediately applicable to assigncontact phases and find positions of assist features. A compositeobjective function is proposed for the non-linear optimizations thatcombines objectives of image fidelity, contour fidelity, and penalizednon-smooth and out of tone solutions. We applied method of localvariations and a gradient descent to the non-linear problem. We proposedelectrical field caching technique. Significant speedup is achieved inthe descent algorithms by using analytical gradient of the objectivefunction. This enables layout inversion on large scale as M log Moperation for M pixels.

Still further mathematical detail of a method of calculating mask pixeltransmission characteristics in accordance with an embodiment of thepresent invention is set forth in U.S. Provisional Patent ApplicationNo. 60/657,260, which is incorporated by reference herein as well as isin the paper “Solving Inverse Problems of Optical Microlithography” byYuri Granik of Mentor Graphics Corporation, reproduced below (withslight edits).

The direct problem of microlithography is to simulate printing featureson the wafer under given mask, imaging system, and processcharacteristics. The goal of inverse problems is to find the best maskand/or imaging system and/or process to print the given wafer features.In this study we will describe and compare solutions of inverse maskproblems.

Pixel-based inverse problem of mask optimization (or “layout inversion”)is harder than inverse source problem, especially for partially-coherentsystems. It can be stated as a non-linear constrained minimizationproblem over complex domain, with large number of variables. We comparemethod of Nashold projections, variations of Fienap phase-retrievalalgorithms, coherent approximation with deconvolution, local variations,and descent searches. We propose electrical field caching technique tosubstantially speedup the searching algorithms. We demonstrateapplications of phase-shifted masks, assist features, and masklessprinting.

We confine our study to the inverse problem of finding the best mask.Other inverse problems like non-dense mask optimization or combinedsource/mask optimization, however important, are not scoped. We alsoconcentrate on the dense formulations of problems, where mask isdiscretized into pixels, and mostly skip the traditional edge-based OPC[25] and source optimization approaches [1].

The layout inversion goal appears to be similar or even the same asfound in Optical Proximity Correction (OPC) or Resolution EnhancementTechniques (RET). However, we would like to establish the inverse maskproblem as a mathematical problem being narrowly formulated, thoroughlyformalized, and strictly solvable, thus differentiating it from theengineering techniques to correct (“C” in OPC) or to enhance (“E” inRET) the mask. Narrow formulation helps to focus on the fundamentalproperties of the problem. Thorough formalization gives opportunity tocompare and advance solution techniques. Discussion of solvabilityestablishes existence and uniqueness of solutions, and guidesformulation of stopping criteria and accuracy of the numericalalgorithms.

The results of pixel-based inversions can be realized by the opticalmaskless lithography (OML) [31]. It controls pixels of 30×30 nm (inwafer scale) with 64 gray levels. The mask pixels can also have negativereal values, which enables phase-shifting.

Strict formulations of the inverse problems, relevant to themicrolithography applications, first appear in pioneering studies ofB.E.A. Saleh and his students S. Sayegh and K. Nashold. In [32], Sayeghdifferentiates image restoration from the image design (a.k.a. imagesynthesis). In both, the image is given and the object (mask) has to befound. However, in image restoration, it is guaranteed that the image isachieved by some object. In image design the image may not be achievableby any object, so that we have to target the image as close as possibleto the desired ideal image. The difference is analogical to solving fora function zero (image restoration) and minimizing a function (imagedesign). Sayegh proceeds to state the image design problem as anoptimization problem of minimizing the threshold fidelity error F_(c) intrying to achieve the given threshold θ at the boundary C of the targetimage ([32], p. 86):

$\begin{matrix}{{{F_{C}\lbrack {m( {x,y} )} \rbrack} = {\oint\limits_{C}{( {{I( {x,y} )} - \theta} )^{n}{{l\min}}}}},} & ( {1A} )\end{matrix}$

where n=2 and n=4 options were explored; I(x, y) is image from the maskm(x, y); x, y are image and mask coordinates. Optical distortions weremodeled by the linear system of convolution with a point-spread functionh(x, y), so that

I(x,y)=h(x,y)*m(x,y),  (2A)

and for the binary mask

m(x,y)={0,1}.  (3A)

Sayegh proposes algorithm of one at a time “pixel flipping”. Mask isdiscretized, and then pixel values 0 and 1 are tried. If the error (1)decreases, then the pixel value is accepted, otherwise it is rejected,and we try the next pixel.

Nashold [22] considered a bandlimiting operator in the place of thepoint-spread function (2A). Such formulation facilitates application ofthe alternate projection techniques, widely used in image processing forthe reconstruction and is usually referenced as Gerchberg-Saxton phaseretrieval algorithm [7]. In Nashold formulation, one searches for acomplex valued mask that is bandlimited to the support of thepoint-spread function, and also delivers images that are above thethreshold in the bright areas B and below the threshold in the darkareas D of the target:

x,y εB:I(x,y)>θ,

x,y εD:I(x,y)<θ.  (4A)

Both studies [32] and [22] advanced solution of inverse problems for thelinear optics. However, the partially coherent optics ofmicrolithography is not a linear but a bilinear system [29], so thatinstead of (2A) the following holds:

I(x,y)=∫∫∫∫q(x−x ₁ ,x−x ₂ ,y−y ₁ ,y−y ₂)m(x ₁ ,y ₁)m*(x ₂ ,y ₂)dx ₁ dx ₂dy ₁ dy ₂,  (5A)

where q is a 4D kernels of the system. While the pixel flipping [32] isalso applicable to the bilinear systems, Nashold technique relies on thelinearity. To get around this limitation, Pati and Kailath [25] proposedto approximate bilinear operator by one coherent kernel h, thepossibility that follows from Gamo results [6]:

I(x,y)≈λ|h(x,y)*m(x,y)|²,  (6A)

where constant λ is the largest eigenvalue of q, and h is thecorrespondent eigenfunction. With this the system becomes linear in thecomplex amplitude A of the electrical field

A(x,y)=√{square root over (λ)}h(x,y)*m(x,y).  (7A)

Because of this and because h is bandlimited, the Nashold technique isapplicable.

Y. Liu and A. Zakhor [19, 18] advanced along the lines started by thedirect algorithm [32]. In [19] they introduced optimization objective asa Euclidean distance ∥•∥₂ between the target I_(ideal) and actual waferimages

F _(I) [m(x,y)]=∥I(x,y)−I _(ideal)(x,y)∥₂→min.  (8A)

This was later used in (1A) as image fidelity error in sourceoptimization. In addition to the image fidelity, the study [18]optimized image slopes in the vicinity of the target contour C :

$\begin{matrix}{{{F_{S}\lbrack {m( {x,y} )} \rbrack} = {{\oint\limits_{C - ɛ}{{I( {x,y} )}{l}}} - {\oint\limits_{C + ɛ}{{I( {x,y} )}{{l\min}}}}}},} & ( {9A} )\end{matrix}$

where C+c is a sized up and C−ε is a sized down contour C ; c is a smallbias. This objective has to be combined with the requirement for themask to be a passive optical element m(x,y)m*(x,y)≦1 or, using infinitynorm ∥.∥_(∞)=max|.|, we can express this as

∥m(x,y)∥_(∞)≦1.  (10A)

In case of the incoherent illumination

I(x,y)=h(x,y)²*(m(x,y)m*(x,y))  (12A)

the discrete version of (9A, 10A) is a linear programming (LP) problemfor the square amplitude p_(i)=m_(i)m_(i)* of the mask pixels, and wasaddressed by the “branch and bound” algorithm. When partially coherentoptics (4A) is considered, the problem is complicated by theinteractions m_(i)m_(j)* between pixels and becomes a quadraticprogramming (QP) problem. Liu [18] applied simulated annealing to solveit. Consequently, Liu and Zakhor made important contributions to theunderstanding of the problem. They showed that it belongs to the classof the constrained optimization problems and should be addressed assuch. Reduction to LP is possible; however, the leanest relevant tomicrolithography and rigorous formulation must account for the partialcoherence, so that the problem is intrinsically not simpler than QP. Newsolution methods, more sophisticated than the “pixel flipping”, havealso been introduced.

The first pixel-based pattern optimization software package wasdeveloped by Y.-H. Oh, J-C Lee, and S. Lim [24], and called OPERA, whichstands for “Optical Proximity Effect Reducing Algorithm.” Theoptimization objective is loosely defined as “the difference between theaerial image and the goal image,” so we assume that some variant of (7A)is optimized. The solution method is a random “pixel flipping”, whichwas first tried in [32]. Despite the simplicity of this algorithm, itcan be made adequately efficient if image intensity can be quicklycalculated when one pixel is flipped. The drawback is that pixelflipping can easily get stuck in the local minima, especially for PSMoptimizations. In addition, the resulting patterns often have numerousdisjoined pixels, so they have to be smoothed, or otherwisepost-processed, to be manufacturable [23]. Despite these drawbacks, ithas been experimentally proven in [17] that the resulting masks can bemanufactured and indeed improve image resolution.

The study [28] of Rosenbluth, A., et al., considered mask optimizationas a part of the combined source/mask inverse problem. Rosenbluthindicates important fundamental properties of inverse mask problems,such as non-convexity, which causes multiple local minima The solutionalgorithm is designed to avoid local minima and is presented as anelaborate plan of sequentially solving several intermediate problems.

Inspired by the Rosenbluth paper and based on his dissertation and theSOCS decomposition [3], Socha delineated the interference mappingtechnique [34] to optimize contact hole patterns. The objective is tomaximize sum of the electrical fields A in the centers (x_(k), y_(k)) ofthe contacts k=1 . . . N :

$\begin{matrix}{{F_{B}\lbrack {m( {x,y} )} \rbrack} =  {- {\sum\limits_{k}{A( {x_{k},y_{k}} )}}}arrow{\min.} } & ( {13A} )\end{matrix}$

Here we have to guess the correct sign for each A(x_(k), y_(k)), becausethe beneficial amplitude is either a large positive or a large negativenumber ([34] uses all positive numbers, so that the larger A thebetter). When kernel h of (7A) is real (which is true for theunaberrated clear pupil), A and F_(B) are also real-valued underapproximation (7A) and for the real mask m. By substituting (7A) into(13A), we get

$\begin{matrix}\begin{matrix}{{- {\sum\limits_{k}{A( {x_{k},y_{k}} )}}} =  {- {\sum\limits_{k}( {h*m} )}} |_{{x = x_{k}},{y = y_{k}}}} \\{= {{\sum\limits_{k}{( {h*m} ) \cdot {\delta ( {{x - x_{k}},{x - x_{k}}} )}}} =}} \\{{= {{- ( {h*m} )} \cdot ( {\sum\limits_{k}{\delta ( {{x - x_{k}},{x - x_{k}}} )}} )}},}\end{matrix} & ( {14A} )\end{matrix}$

where the dot denotes an inner product ƒ·g=∫∫fgdxdy. Using the followingrelationship between the inner product, convolution* , andcross-correlation ∘ of real functions

(ƒ*g)·p=ƒ·(g∘p),  (15A)

$\begin{matrix}{{{- {\sum\limits_{k}{A( {x_{k},y_{k}} )}}} = {{{- ( {h \circ {\sum\limits_{k}{\delta ( {{x - x_{k}},{x - x_{k}}} )}}} )} \cdot m} = {{- G_{b}} \cdot m}}},} & ( {16A} )\end{matrix}$

where function G_(I) is the interference map [34]. With (16A) theproblem (13A) can be treated as LP with simple bounds (as defined in[8]) for the mask pixel vector m={m_(i)}

−G _(b) ·m→min

−1≦m_(i)≦1.  (17A)

In an innovative approach to the joined mask/source optimization byErdmann, A., et al. [4], the authors apply genetic algorithms (GA) tooptimize rectangular mask features and parametric source representation.GA can cope with complex non-linear objectives and multiple localminima. It has to be proven though, as for any stochastically basedtechnique that the runtime is acceptable and quality of the solutions isadequate. Here we limit ourselves to the dense formulations and moretraditional mathematical methods, so the research direction of [4] and[15] however intriguing is not pertinent to this study.

The first systematic treatment of source optimization appeared in [16].This was limited to the radially-dependent sources and periodic maskstructures, with the Michelson contrast as an optimization objective.Simulated annealing is applied to solve the problem. After this study,parametric [37], contour-based [1], and dense formulations [28], [12],[15] were introduced. In [12], the optimization is reduced to solving anon-negative least square (NNLS) problem, which belongs to the class ofthe constrained QP problems. The GA optimization was implemented in [15]for the pixelized source, with the objective to maximize image slopes atthe important layout cutlines.

Reduction to Linear Problem

The inverse mask problem can be reduced to a linear problem, includingtraditional LP, using several simplification steps. The first step is toaccept coherent approximation (6A,7A). Second, we have to guesscorrectly the best complex amplitude A_(ideal) of the electrical fieldfrom

I_(ideal)=A_(ideal)A_(ideal)*  (18A)

where I_(ideal) is the target image. If we consider only the real masksm=Re[m] and real kernels h=Re[h], then from (7A) we conclude that A isreal and thus we can set A_(ideal) to be real-valued. From (18A) we get

A _(ideal) =±√{square root over (I_(ideal))},  (19A)

which means that A_(ideal) is either +1 or −1 in bright areas of thetarget, and 0 in dark areas. If the ideal image has M bright pixels, thenumber of possible “pixel phase assignments” is exponentially large2^(M). This can lead to the phase-edges in wrong places, but of coursecan be avoided by assigning the same value to all pixels within a brightfeature: for N bright features we get 2^(N) different guesses. After wechoose one of these combinations and substitute it as A_(ideal) into(7A), we have to solve

A _(ideal)(x,y)=√{square root over (λ)}h(x,y)*m(x,y)  (20A)

form. This is a deconvolution problem. Within the zoo of deconvolutionalgorithms, we demonstrate Weiner filtering, which solves (20A) in someleast square sense. After applying Fourier transformation F[ . . . ] to(20A) and using convolution theorem F[h*m]=F[h]F[m], we get

$\begin{matrix}{{\hat{m} = {\frac{{\hat{A}}_{ideal}}{\sqrt{\lambda}\hat{h}} = {{\hat{A}}_{ideal}\frac{{\hat{h}}^{*}}{\sqrt{\lambda}{\hat{h}}^{2}}}}},} & ( {21A} )\end{matrix}$

where the circumflex denotes Fourier transforms: {circumflex over(m)}=F[m], ĥ=F[h]. The Wiener filter is a modification of (21A) where arelative noise power P is added to the denominator, which helps to avoiddivision by 0 and suppresses high harmonics:

$\begin{matrix}{\hat{m} = {{\hat{A}}_{ideal}{\frac{{\hat{h}}^{*}}{{\sqrt{\lambda}{\hat{h}}^{2}} + P}.}}} & ( {22A} )\end{matrix}$

Final mask is found by the inverse Fourier transformation:

$\begin{matrix}{m = {{F^{- 1}\lbrack {{\hat{A}}_{ideal}\frac{{\hat{h}}^{*}}{{\sqrt{\lambda}{\hat{h}}^{2}} + P}} \rbrack}.}} & ( {23A} )\end{matrix}$

As the simplest choice we set P=const>0 to be large enough to satisfymask constraint (11A). The results are presented in FIGS. 11A-11D. Theimaging is simulated for annular 0.7/0.5 optical model with 0.75 NA and193 nm lambda. The first inset on the left shows contours when mask isthe same as target. The space between horizontal lines and L-shape tendsto bridge, and the comb structure is pinching. The contour fidelityafter deconvolution is much better (right inset)., the bridging andpinching tendencies are gone. The semi-isolated line width is on target.The contrast is also improved, especially for the semi-isolated verticalline. However, the line ends and corner fidelity is not improved. Themask contours after deconvolution are shown in the right bottom inset.Positive transmission is assigned to all main features, which arecanvassed with the assist features of negative transmission.

We can also directly solve (2) in the least square sense

∥√{square root over (λ)}h(x,y)*m(x,y)−A _(ideal)(x,y)∥→min.  (24A)

In the matrix form

∥Hm−A _(ideal)∥→min

H_(ij)=√{square root over (λ)}h_(i−j).  (25A)

Matrix H has multiple small eigenvalues. The problem is ill-posed. Thestandard technique dealing with this is to regularize it by adding normof the solution to the minimization objective [14]:

∥Hm −A _(ideal)∥² +α∥m∥ ²→min,  (26A)

where the regularization parameter α is chosen from secondaryconsiderations. In our case we chose α large enough to achieve∥m∥_(∞)=1. The problem (26A) belongs to the class of unconstrainedconvex quadratic optimization problems, with guaranteed unique solutionin non-degenerate cases. It can be solved by the methods of linearalgebra, because (26) is equivalent to solving

(H+αI)m=A _(ideal)  (27A)

by the generalized inversion [12] of the matrix H+αI. The results arepresented in FIGS. 12A-12D.

This method delivers pagoda-like corrections to image corners. Somehints of hammer-heads and serifs can be seen in mask contours. Line endsare not corrected. Comparison of contrasts between the case when mask isthe same as target show improved contrast, especially between the comband semi-isolated line.

Further detailing of the problem (26A) is possible by explicitly addingmask constrains, that is we solve

∥Hm−A _(ideal)∥² −α∥m∥ ²→min

∥m∥_(∞)≦1,  (28A)

This is a constrained quadratic optimization problem. It is convex asany linear least square problem with simple bounds. Convexity guaranteesthat any local minimizer is global, so that the solution method is notconsequential: all proper solvers converge to the same (global)solution. We used MATLAB routine Isqlin to solve (29A). The results arepresented in FIGS. 13A-13D. This solution has better contrast than thetwo previous ones, with more complex structure of the assist features.

Any linear functional of A is also linear by m, in particular we canform a linear objective by integrating A over some parts of the layout,as in (13A). One of the reasonable objectives to be formed by suchprocedure is the sum of electrical amplitudes through the region B,which consists of the all or some parts of the bright areas:

$\begin{matrix}{{{F_{B}\lbrack {m( {x,y} )} \rbrack} =  {- {\int{\int_{B}{{A( {x,y} )}\ {x}{y}}}}}arrow\min },} & ( {28\; B} )\end{matrix}$

that is we try to make bright areas as bright as possible. Using thesame mathematical trick as in (14A), this is reduced to the linearobjective

−G _(b) ·m→min,  (29A)

where G_(b)=h∘b, and b is a characteristic function of the bright areas.This seems to work well as the basis for the contact optimizations. Itis harder to form region B for other layers. If we follow suggestion [4]to use centers of the lines, then light through the corners becomesdominant, spills over to the dark areas, and damages image fidelity.This suggests that we have to keep dark areas under control as well.Using constraints similar to (4A), we can require for each dark pixel tobe of the limited brightness θ

x,y εD:−θ≦A(x,y)≦θ,  (30A)

or in the discrete form

−θ≦H_(d)m≦θ.  (30B)

where H_(d) is matrix H without rows correspondent to the brightregions. Though equations (28A) and (30B) form a typical constrained LPproblem, MATLAB simplex and interior point algorithms failed toconverge, perhaps because the matrix of constraints has largenull-space.

The linearization (7A) can be augmented by the threshold operator tomodel the resist response. This leads to Nashold projections [22].Nashold projections belong to the class of the image restorationtechniques, rather than to the image optimizations, meaning that themethod might not find the solution (because it does not exists at all),or in the case when it does converge, we cannot state that this solutionis the best possible. It has been noted in [20] that the solutionsdepend on the initial guess and do not deliver the best phase assignmentunless the algorithm is steered to it by a good initial guess. Moreover,if the initial guess has all phases set to 0, then so has the solution.

Nashold projections are based on Gerchberg-Saxton [7] phase retrievalalgorithm. It updates a current mask iterate m^(k) via

m ^(k+1)=(P _(m) P _(s))m ^(k),  (31A)

where P_(s) is a projection operator into the frequency support of thekernel h, and P_(m) is a projection operator that forces thethresholding (4A). Gerchberg-Saxton iterations tend to stagnate. Fienap[5] proposed basic input-output (BIO) and hybrid input-output (HIO)variations that are less likely to be stuck in the local minima Thesevariations can be generalized in the expression

m ^(k+1)=(P _(m) P _(s)+α(γ(P _(m) P _(s) −P _(s))−P _(m) +I)m^(k),  (32A)

where I is an identity operator; α=1,γ=0 for BIO, α=1,γ=1 for HIO, andα=0,γ=0 for the Gerchberg-Saxton algorithm.

We implemented operator P. as a projection onto the ideal image

$\begin{matrix}{{{P_{m}m^{k}} = {\frac{m^{k}}{m^{k}}\sqrt{I_{ideal}}}},} & ( {33A} )\end{matrix}$

and P_(s) as a projection to the domain of the kernel h, i.e. P_(s)zeros out all frequencies of {circumflex over (m)} which are high thanthe frequencies of the kernel h. The iterates (32A) are very sensitiveto the values of its parameters and the shape of ideal image. We wereable to find solutions only when the ideal image is smoothed. We usedGaussian kernel with the diffusion length of 28 nm, which is slightlylarger than the pixel size 20 nm in our examples. The behavior ofiterates (32A) is not yet sufficiently understood [36], whichcomplicates choice of α, γ. We found that in our examples theconvergence is achieved for α=0.9,γ=1 after 5000 iterations. Whenα=0,γ=0, which corresponds to (31), the iterations quickly stagnateconverging to a non-printable mask.

As shown in FIG. 1B, Nashold projections assign alternating phases tothe main features and insert assists between lines. The lines widths areon-target, but line-ends are not corrected. The solution has very goodcontrast. When projections stagnate, the phases alternate along thelines. This “phase entanglement” can sometimes happen in the non-linearproblems (considered in a section below) when their iterations startfrom the random pixel assignments.

Quadratic Problems

In the quadratic formulations of the inverse problems, the coherentlinearization (6A) is not necessarily. We can directly use bilinearintegral (5A). Our goal here is to construct an objective function as isa quadratic form of mask pixels. We start with (8A) and replaceEuclidean norm (norm 2) with Manhattan norm (norm 1):

F _(I) [m(x,y)]=∥I(x,y)−I _(ideal)(x,y)∥₁→min  (34A)

The next step is to assume that the ideal image is sharp, 0 in darkregions and 1 in bright regions, so that I(x,y)≧I_(ideal)(x,y) in thedark regions and I(x,y)≦I_(ideal)(x,y) in the bright regions. This letsus to remove the module operation from the integral (34A):

∥I(x,y)−I _(ideal)(x,y)∥₁ =∫∫|I−I _(ideal) |dxdy=∫∫w(x,y)(I(x,y)−I_(ideal)(x,y))dxdy,  (35A)

where w(x, y) is 1 in dark regions and −1 in bright regions. Finally wecan ignore the constant term in (35A), which leads to the objective

F _(w[m)(x,y)]=∫∫wI(x,y)→min.  (36A)

The weighting function w can be generalized to have any positive valuein dark regions, any negative value in bright regions, and 0 in theregions which we choose to ignore. Proper choice of this function coversthe image slope objective (9A), but not the threshold objective (1A).Informally speaking, we seek to make bright regions as bright aspossible, and dark regions as dark as possible. Substituting (5A) into(36A), we get

∫∫wI(x,y)dxdy =∫∫∫∫Q(x ₁ ,y ₁,x₂ y ₂)m(x ₁ ,y ₁)m*(x ₂ ,y ₂)dx ₁ dx ₂ dy₁ dy ₂,  (37A)

Where

Q(x ₁ ,y ₁ ,x ₂ y ₂)=∫∫w(x,y)q(x−x ₁ ,x−x ₂ ,y−y ₁ ,y−y ₂)dxdy.  (38A)

Discretization of (37A) results to the following constrained QP

F _(w) [m]=m*Qm→min

∥m∥_(∞)≦1.  (39A)

The complexity of this problem depends on the eigenvalues of matrix Q.When all eigenvalues are non-negative, then it is convex QP and anylocal minimizer is global. This is a very nice property, because we canuse any of the numerous QP algorithms to find the global solution and donot have to worry about local minima Moreover, it is well known that theconvex QP can be solved in polynomial time. The next case is when alleigenvalues are non-positive, a concave QP. If we remove constraints,the problem becomes unbounded (no solutions). This means that theconstraints play a decisive role: all solutions, either local or global,end up at some vertex of the box ∥m∥_(∞)≦1. In the worst case scenario,the solver has to visit all vertices to find the global solution, whichmeans that the problem is NP-complete, i.e., it may take an exponentialamount of time to arrive at the global minima The last case is anindefinite QP when both positive and negative eigenvalues are present.This is the most complex and the most intractable case. An indefinite QPcan have multiple minima, all lie on the boundary.

We conjecture that the problem (39A) belongs to the class of indefiniteQP. Consider the case of the ideal coherent imaging, when Q is adiagonal matrix. Vector w lies along its diagonal. This means thateigenvalues μ₁, μ₂ . . . of Q are the same as components of the vectorw, which are positive for dark pixels and negative for bright pixels. Ifthere is at least one dark and one bright pixel, the problem isindefinite. Another consideration is that if we assume that (39A) isconvex, then the stationary internal point m=0 where the gradient iszero

$\begin{matrix}{\frac{{F_{w}\lbrack m\rbrack}}{m} = {{2Q\; m} = 0}} & ( {40A} )\end{matrix}$

is the only solution, which is a trivial case of mask being dark. Thismeans that (39) is either has trivial (global) solution, or it isnon-convex.

Related to (39) QP was considered by Rosenbluth [28]:

m*Q _(d) m→min

m*Q _(b) m→b,  (41A)

where Q_(d) and Q_(b) deliver average intensities in bright and darkregions correspondingly. The objective is to keep dark regions as darkas possible while maintaining average intensity not worse than somevalue b in bright areas. Though the problem was stated for the specialcase of the off-centered point-source, the structure of (41A) is verysimilar to (39A). Using Lagrange multipliers, we can convert (41A) to

m*(Q _(d) −λQ _(b))m→min

∥m∥_(∞)≦1

λ≧0,  (42A)

which is similar to (39A).

Another metric of the complexity of (39A) is number of the variables,i.e., the pixels in the area of interest. According to Gould [10], theproblems with order of 100 variables are small, more than 10³ are large,and more than 10⁵ are huge. Considering that the maskless lithographycan control transmission of the 30 nm by 30 nm pixel [31], the QP (39A)is large for the areas larger than 1 um by 1 um, and is huge for theareas lager than 10 um by 10 um. This has important implications for thetype of the applicable numerical methods: in large problems we can usefactorizations of matrix Q, in huge problems factorizations areunrealistic.

For the large problems, when factorization is still feasible, a dramaticsimplification is possible by replacing the infinity norm by theEuclidean norm in the constraint of (39A), which results in

F_(w) [m]=m*Qm→min

∥m∥₂≦1  (43A)

Here we search for the minimum inside a hyper-sphere versus a hyper-cubein (39A). This seemingly minor fix carries the problem out of the classof NP-complete to P (the class of problems that can be solved inpolynomial time). It has been shown in [35] that we can find globalminima of (43A) using linear algebra. This result served as a base forthe computational algorithm [1] which specifically addresses indefiniteQP.

The problem (43A) has the following physical meaning: we optimize thebalance of making bright regions as bright as possible and dark regionsas dark as possible while limiting light energy ∥m∥hd 2 ² coming throughthe mask. To solve this problem, we use procedures outlined in [31,32].First we form Lagrangian function of (43A)

L(m,λ)=m*Qm+λ(∥m∥ ²−1).  (44A)

From here we deduce the first order necessary optimality conditions ofKarush-Kuhn-Tucker (or KKT conditions, [20]):

2(Q+λI)m=0

λ(∥m∥−1)=0

λ≧0

∥m∥≦1.  (45A)

Using Sorensen [35], we can state what that (43A) has a global solutionif we can find such λ and m that (45A) is satisfied and the matrix Q+λIis positive semidefinite or positively defined. Let us find thissolution.

First we notice that we have to choose λ large enough to compensate thesmallest (negative) eigenvalue of Q, i.e.

λ≧|μ₁|>0.  (46A)

From the second condition in (45A) we conclude ∥m∥=1, that is thesolution lies on the surface of hyper-sphere and not inside it. The lastequation to be satisfied is the first one from (45A). It has anon-trivial ∥m∥>0 solution only when the lagrange multiplier λ equals toa negative of one of the eigenvalues λ=−μ_(i). This condition and (46A)has a unique solution λ=−μ₁, because other eigenvalues μ₂, μ₃, . . . areeither positive so that λ≧0 does not hold, or they are negative, butwith absolute value that is smaller than μ₁, so that does not hold.

After we determined that λ=−μ₁, we can find m from 2(Q−μ₁I)m=0 as thecorresponding eigenvector m=v₁. This automatically satisfies ∥m∥=1,because all eigenvectors are normalized to have a unit length. Weconclude that (43A) has a global solution which corresponds to thesmallest negative eigenvalue of Q. This solution is a good candidate fora starting point in solving (39A): we start from the surface of thehyper-sphere and proceed with some local minimization technique to thesurface of the hyper-cube.

As we have shown, the minimum eigenvalue of Q and its eigenvector playspecial role in the problem by defining the global minimum. However,other negative eigenvectors are also important, because it is easy tosee that any pair

λ=−μ_(i)>0

m=v_(i)  (47A)

is a KKT point and as such defines a local minimum. The problem has asmany local minima as negative eigenvalues. We may also consider startingour numerical minimization from one of these “good” minima, because itis possible that a local minimum leads to a better solution in thehyper-cube than a global minimum of the spherical problem.

FIGS. 2A, 2B and 2C show three strongest local minima of the problem(39A) for the structure of FIG. 1B. These local minima point to thebeneficial interactions between layout features, suggesting alternatingphase assignments. For example, the second solution suggests that theL-shape transmission should be chosen positive, while the comb hasnegative transmission, the dense vertical line of the comb has positivetransmission, and the second horizontal line has negative transmission.

Results of the similar analysis for the case of the contact holes aredisplayed in FIGS. 3B, 3C and 3D. These results are much stronger, andcan be used directly in applications. The method “proposes” beneficialphases for the contacts and position and phases of the assists. The mostinteresting solution is shown in the low right inset, where all contactshave well-defined transmissions, with 3 contacts positive and 4 contactsnegative. The advantages of this method comparing to with IML [3] isthat this method automatically finds the best phases of the contacts andis not based of the coherent approximation.

FIGS. 3A-3D shows the first three local minima for QP on a hyper-spherefor the contact holes and process conditions from Socha [3]. The thirdsolution has the clearest phase assignments and the position of assists.

For the positive masks, in particular for the binary masks, theconstraint can be tightened to ∥m−0.5∥_(∞)≦0.5. Then the correspondentto (39A) problem is

F _(w) [m]=m*Qm→min

∥m−0.5∥_(∞)≦0.5  .(48A)

This is also an indefinite QP and is NP-complete. Replacing hereinfinity norm with Euclidean norm, we get a simpler problem

m*Qm→min

∥Δm∥₂≦0.5

Δm=m−m ₀ ,m ₀={0.5,0.5, . . . , 0.5}.  (49A)

The Lagrangian can be written as

L(m,λ)=m*Qm+λ(∥m−m ₀∥²−0.25).  (50A)

The KKT point must be found from the following conditions

(Q+λI)Δm=−Qm ₀

λ(∥Δm∥ ²−0.25)=0

λ≧0

∥m∥≦0.5.  (51A)

This is more complex problem than (45A) because the first equation isnot homogeneous and the pairs λ=−μ_(i),Δm=v_(i) are clearly not thesolutions. We can still apply the condition of the global minimumλ≧−μ₁>0 (Sorensen [35]). From the second condition we conclude that∥Δm∥²=0.25, meaning that all solutions lie on the hyper-sphere with thecenter at m₀. The case λ=−μ₁ is eliminated because the first equation isnot homogeneous, so that we have to consider only λ>−μ₁. Then Q+λI isnon-singular, we can invert it, and find the solution

Δm=−(Q+λI)⁻¹ Qm ₀.  (52A)

The last step is to find the lagrange multiplier A that satisfy theconstraint ∥Δ∥²=0.25, that is we have to solve

∥(Q+λI)⁻¹ Qm ₀∥=0.5  .(53A)

This norm monotonically increases from 0 to infinity in the interval−∞<λ<−μ₁, thus (53A) has to have exactly one solution in this interval.The pair A, Am that solves (52A-53A) is a global solution of (49A). Weconjecture that there are fewer KKT points of local minima of (49A) thanin (45A) (may be there are none), but this remains to be proven byanalyzing behavior of the norm (53A) when lagrange multiplier is betweennegative eigenvalues. The solutions of (49A) show how to insert assistfeatures when all contacts have the same phases.

General Non-Linear Problems

Consider objective (8A) of image fidelity error

F _(I) [m(x,y)]=∥(x,y)−I _(ideal)(x,y)∥→min.  (54A)

We can state this in different norms, Manhattan, infinity, Euclidean,etc. The simplest case is a Euclidean norm, because (54A) becomes apolynomial of the forth degree (quartic polynomial) of mask pixels. Theobjective function is very smooth in this case, which ease applicationof the gradient-descent methods. While theory of QP is well developed,the polynomial optimization is an area of growing research interest, inparticular for quartic polynomials [27].

We can generalize (54A) by introducing weighting w=w(x, y) to emphasizeimportant layout targets and consider smoothing in Sobolev norms as in[12]:

F _(w) [m(x,y)]²=∥√{square root over (w)}·(I−I _(ideal))₂ ²+α₁ ∥L ₁ m∥ ₂²+α₃ ∥m−m ₀∥₂ ²→min,  (55A)

where L₁,L₂ are the operators of first and second derivatives, m₀=m₀(x,y) is some preferred mask configuration that we want to be close to (forexample, the target), and α₁, α₂, α₃ are smoothing weights. Thesolutions of (55A) increase image fidelity; however, the numericalexperiments show that the contour fidelity of the images is notadequate. To address, we explicitly add (1A) into (55A):

$\begin{matrix}{{F_{wc}\lbrack m\rbrack}^{2} =  {{{\sqrt{w} \cdot ( {I - I_{ideal}} )}}_{2}^{2} + {\oint\limits_{C}{( {I - \theta} )^{n}{l}}} + {\alpha_{1}{{L_{1}m}}_{2}^{2}} + {\alpha_{2}{{L_{2}m}}_{2}^{2}} + {\alpha_{3}{{m - m_{0}}}_{2}^{2}}}arrow{\min.} } & ( {56A} )\end{matrix}$

If the desired output is a two-, tri-, any other multi-level tone mask,we can add penalty for the masks with wrong transmissions. The simplestform of the penalty is a polynomial expression, so for example for thetri-tone Levenson-type masks with transmissions −1, 0, and 1, weconstruct the objective as

$\begin{matrix}{{{F_{wce}\lbrack m\rbrack}^{2} =  {{{{\sqrt{w} \cdot ( {I - I_{ideal}} )}}_{2}^{2} + {\oint\limits_{C}{( {I - \theta} )^{n}{{l++}}\alpha_{1}{{L_{1}m}}_{2}^{2}}} + {\alpha_{2}{{L_{2}m}}_{2}^{2}} + {\alpha_{3}{{{m - m_{0}}}_{2}^{2}++}\alpha_{4}{{( {m + e} ){m( {m - e} )}}}_{2}^{2}{m}_{\infty}}} \leq 1}arrow\min },} & ( {57A} )\end{matrix}$

where e is a one-vector. Despite all the complications, the objectivefunction is still a polynomial of the mask pixels. To optimize for thefocus depth, the optimization of (57A) can be conducted off-focus, aswas suggested in [16, 12]. After discretization, (55A) becomes anon-linear programming problem with simple bounds.

We expect that this problem inherits property of having multiple minimafrom the corresponding simpler QP, though smoothing operators of (57A)have to increase convexity of the objective. In the presence of multiplelocal minima the solution method and staring point are highlyconsequential: some solvers tend to converge to the “bad” localsolutions with disjoined masks pixels and entangled phases, othersbetter navigate solution space and chose smoother local minima TheNewton-type algorithms, which rely on the information about secondderivatives, should be used with a caution, because in the presence ofconcavity in (57A), the Newtonian direction may not be a descentdirection. The branch-and-bound global search techniques [18] are notthe right choice because they are not well-suited for the largemulti-dimensional optimization problems. Application of stochastictechniques of simulated annealing [24] or GA [4] seems to be anoverkill, because the objective is smooth. It is also tempting toperform non-linear transformation of the variables to get rid of theconstraints and convert problem to the unconstrained case, for exampleby using transformation x_(i)=tanh(m_(i)) or m_(i)=sin(x_(i)), however,this generally is not recommended by experts [8, p. 267].

The reasonable choices to solve (57A) are descent algorithms withstarting points found from the analytical solutions of the related QP.We apply algorithms of local variations (“one variable at a time”),which is similar in spirit to the pixel flipping [32, 24], and also usea variation of the steepest descent by Frank and Wolfe [21] to solveconstrained optimization problems.

In the method of local variation, we chose the step Δ₁ to compare threeexploratory transmissions for the pixel i: m_(i) ¹, m_(i) ¹+Δ₁, andm_(i) ¹−Δ₁. If one of these values violates constraints, then it ispulled back to the boundary. The best of these three values is accepted.We try all pixels, optionally in random exhaustive or circular order,until no further improvement is possible. Then we reduce step Δ₂<Δ₁ andrepeat the process until the step is deemed sufficiently small. Thisalgorithm is simple to implement. It naturally takes care of the simple(box) constraints and avoids the general problem of other moresophisticated techniques (like Newton), which may converge prematurelyto a non-stationary point. This algorithm calculates the objectivefunction numerous times; however, the runtime cost of its exploratorycalls is very low with the electrical field caching (see the nextsection). Other algorithms require fewer but more costly non-exploratorycalls. This makes method of local variation a legitimate tool in solvingthe problem.

Frank and Wolfe method is an iterative algorithm to solve constrainproblems. At each step k we calculate the gradient ∇F^(k) of theobjective and then replace the non-linear objective with its linearapproximation. This reduces the problem to LP with simple bounds:

∇F ^(k) ·m→min

∥m∥_(∞)≦1.  (59A)

The solution of this m=l^(k) is used to determine the descent direction

p ^(k)=l^(k) −m ^(k−1).  (60B)

Then the line search is performed in the direction of p^(k) to minimizethe objective as a function one variable γε[0,1]:

F[m ^(k−1) +γp ^(k)]→min.  (61B)

The solution m^(k)=m^(k−1)+γp^(k)is accepted as the next iterate. Theiterations continue until convergence criteria are met. Electrical fieldcaching helps to speedup the gradient calculations and line search ofthis procedure.

In FIGS. 14A-14C we show results of solving (55A) for the positive mask0≦m1. The assist features can be seen around main structures. However,the solution comes up very bright and the contrast is only marginallyimproved. This is corrected in (56A) with the introduction of thecontour fidelity metric. The results are shown in FIGS. 15A-15C.

FIGS. 15A-15C illustrates a method of local variations for the objective(57A) with contour fidelity. The contours are on target almosteverywhere, including line ends. The image contrast is improved. Maskhas rows of assist features and serifs.

Result of the local variation algorithm for the PSM mask are shown inFIG. 16. The contours are on target, except some corner rounding. Themain features have the same phase. Each side of the pads and lines isprotected by the queue of assist features with alternating phases.

Next example demonstrate solutions when main features have the samephase and assist features can have phase shift, FIGS. 17A-17C. Weobserve negative transmission of the assists on the mask. The contrastalong the cutline is improved in comparison to the ideal case (maskequal target). Contour fidelity is very good, the third inset. The lastexample is contact holes, FIG. 18. The method is capable of insertingassist contacts and deliver complex interferometric assist features inPSM case.

FIGS. 18A-18F illustrates a contact holes example. First row showsbinary mask. Small assist contact holes are inserted around maincontacts. The image contrast is compared to the case when mask is thesame as target. The contrast is improved significantly. Image contoursare on target, the third column. Second row demonstrates PSM mask, withcomplex structure of assists holes, which are hard to separate from themain features. The contrast is even better than for the binary mask.Despite very complex mask, the contours are on target (lower rightinset) and sidelobs do not print.

Electrical Field Caching

The speed of the descent and local variation algorithms criticallydepends on the ability to quickly re-calculate image intensity when oneor a few pixels change. We use electrical field caching procedure tospeedup this process.

According to SOCS approximation [3], the image intensity is thefollowing sum of convolutions of kernels h_(i)(x, y) with the mask m(x,y):

$\begin{matrix}{{{I( {x,y} )} = {\sum\limits_{i = 1}^{N}{\lambda_{i}{A_{i}( {x,y} )}{A_{i}^{*}( {x,y} )}}}},{A_{i} = {{h_{i}( {x,y} )}*{{m( {x,y} )}.}}}} & ( {60C} )\end{matrix}$

Suppose that we know the electrical fields A_(i) ⁰ for the mask m⁰ andwant to calculate intensity for the slightly different mask m′. Then

A _(i) ′=A _(i) ⁰ +h _(i)*(m′−m ⁰).  (61C)

These convolutions can be quickly calculated by the directmultiplication, which is O(d.M.N) operation, where d is the number ofdifferent pixels between m⁰ and m′, M is pixel count of the kernels, andN is number of kernels. This is faster than convolution by FFT when O(d)is smaller than O(log(M)). Constantly updating the cache A_(i) ⁰, we canquickly re-calculate intensities for small evolutionary mask changes.Formula (61A) is helpful in gradient calculations, because they alterone pixel at a time.

The additivity of the electrical fields can also be exploited to speedupintensity calculations in the line search (61A). If the mask m^(k−1)delivers electrical fields A_(i) ^(k−1), and the mask p^(k) deliversB_(i) ^(k), then the intensity from the mask m=m^(k−1)+γp^(k) can bequickly calculated through its electrical fields A_(i):

A _(i) =A _(i) ^(k−1) +γB _(i) ^(k).  (62A)

This avoids expensive convolutions of (60C).

In one embodiment of the invention, the optimization function to beminimized in order to define the optimal mask data has the form

$\begin{matrix}{{Min}{\sum\limits_{i}{w_{i} \cdot {{I_{i} - I_{i_{1}{ideal}}}}^{2}}}} & ( {64A} )\end{matrix}$

where I_(i)=image intensity evaluated at a location on the wafercorresponding to a particular pixel in the mask data. Typically I_(i)has a value ranging between 1 (full illumination) and φ (noillumination); I_(i) _(l) _(ideal) the desired image intensity on thewafer at a point corresponding to the pixel and w_(i)=a weighting factorfor each pixel.

As indicated above, the ideal image intensity level for any point on thewafer is typically a binary value having an intensity level of 0 forareas where no light is desired on a wafer and 1 where a maximum amountof light is desired on the wafer.

However, in some embodiments, a better result can be achieved if themaximum ideal image intensity is set to a value that is determined fromexperimental results, such as a test pattern for various pitchesproduced on a wafer with the same photolithographic processing systemthat will be used to produce the desired layout pattern in question. Abetter result can also be achieved if the maximum ideal image intensityis set to a value determined by running an image simulation of a testpattern for various pitches, which predicts intensity values that willbe produced on a wafer using the photolithographic processing systemthat will be used in production for the layout pattern in question.

FIG. 21 illustrates a one example of a test pattern for various pitches150 having a number of vertically aligned bars 152, 154, 156, . . . 166.As will be appreciated by those skilled in the art, photolithographicprocessing systems can be characterized by the tightest pattern offeatures that the system can reliably produce on a wafer. The testpattern for various pitches 150 includes features at this tightest pitchand features that are spaced further apart. In one embodiment of theinvention, the maximum ideal image intensity I_(ideal) defined for apoint on a wafer is determined by simulating exposure in thephotolithographic system using the test pattern 150 and determining themaximum image intensity in the area of the tightest pitch features Theminimum ideal image intensity is generally selected to be below theprint threshold of the resist materials to be used. Typically, theminimum ideal image intensity is set to zero.

With the maximum ideal image intensity determined using the test pattern150, the transmission values of the pixels in the mask data that willresult in the objective function (64A) being minimized are thendetermined Once the transmission values of the pixels has beendetermined, the mask pixel data is converted to a suitable mask writingformat, and provided to a mask writer that produces one or masks. Insome embodiments, a desired layout is broken into one or more frames andmask pixel data is determined for each of the frames.

FIG. 22 illustrates one example of an optimized mask pattern that willreproduce a desired pattern of features such as the test pattern 150shown in FIG. 21. The optimized mask data 180 includes a number oflarger mask features 182, 184, 186, . . . 196 that, when exposed, willcreate each of the vertical bars 152-166 in the test pattern. Inaddition, the optimized data 180 includes a number of additionalfeatures 200, 202, etc., surrounding the larger mask features 182-196.The additional features 200, 202 are a result of the mask pixeloptimization process described above. In one embodiment of theinvention, the additional features 200, 202 are simulated on a mask assubresolution assist features (SRAFs) such that they themselves will notform an image that prints in resist on a wafer when exposed duringphotolithographic processing.

FIG. 23 illustrates the simulated image intensity on a wafer whenexposed through a conventional mask using the test pattern 150. Inaddition, FIG. 23 shows a simulation of the image intensity on a waferwhen exposed through a mask having optimized mask pattern data such asshown in FIG. 22. In FIG. 23, the graph includes a number of minima 222,224, 226, . . . 236 where the image intensity is minimized correspondingto each of the vertical bars 152, 154, 156, . . . 166 in the testpattern 150 shown in FIG. 21. Outside the area of the features, theimage intensity increases where more light will reach the wafer. With aconventional mask, the variations in intensity are greater and varybetween approximately 0.05 and 0.6. However, when a wafer is exposedthrough a mask having the optimized mask data such as illustrated inFIG. 22, the variations in intensity are smaller such that imageintensity varies between approximately 0.05 and approximately 0.25. Ascan be seen in FIG. 23, the maximum image intensity obtained from anoptimized mask pattern such as that shown in FIG. 22 has a moreconsistent or uniform image intensity because the mask pattern mimicsthe closely spaced features from the test pattern.

As indicated above, each pixel in the objective function may be weightedby the weight function w. In some embodiments, the weight function isset to 1 for all pixels and no weighting is performed. In otherembodiments, the weights of the pixels can be varied by, for example,allowing all pixels within a predetermined distance from the edge of afeature in the ideal image to be weighted more than pixels that are farfrom the edges of features. In another embodiment, pixels within apredetermined distance from designated tagged features (e.g. featurestagged to be recognized as gates, or as contact holes, or as line ends,etc.) are given a different weight. Weighting functions can be set tovarious levels, depending on the results intended. Typically, pixelsnear the edge of the ideal image would have a weight w=10, while thosefurther away would have a weight w=1. Likewise, line ends (whose imagesarea known to be difficult to form accurately) may be given a smallerweight w=0.1, while other pixels in the image may be given a weight w=1.Both functions may also be applied (i.e. regions near line ends have aweight w=0.1, and the rest of the image has a weight w=1 except near theedges of the ideal image away from the line ends, where the weight wouldbe w=10.) Alternatively, if solutions using SRAFs are desired, and theseSRAFs occur at a predetermined spacing relative to main features,weighting functions which have larger values at locations correspondingto these SRAF positions can be constructed.

It should be noted that the absolute values of weighting functions canbe set to any value; it is the relative values of the weightingfunctions across pixels that makes them effective. Typically, distinctregions have relative values that differ by a factor of 10 or more toemphasize the different weights.

In some instances it has been found that by setting the maximum idealimage intensity for use in an objective function to the maximum imageintensity for tightly spaced features of the test pattern 150, theprocess window of the photolithographic processing is increased.

FIGS. 24A and 24B illustrate two possible techniques for generation ofoptimized mask data on one or more photolithographic masks. As indicatedabove, the optimized mask data may produce a set of mask featurescorresponding to desired layout features to be created on a wafer. Inthe example shown in FIG. 24, the mask data includes a feature 250 thatcorresponds to a square via or other small square feature to be createdon a wafer. However, the group of mask pixels defining the feature 250in the optimized mask data has a generally circular shape. Because maskwriters are not generally able to easily produce curved or circularstructures, the mask data for the mask feature 250 can be simulated witha generally square polygon 260. In addition, the mask data includes anadditional feature 252 that is a result of the optimization process. Theadditional feature 252 surrounds the desired feature 250 and has agenerally annular shape. Again, such a curved annular feature would bedifficult to accurately produce using most mask writers. One techniquefor emulating the effect of the annular feature 252 is to use a numberof rectangular polygons 262, 264, 266, 268 positioned in the area of theannular feature 252. Each of the polygons 262-268 has a size that isselected such that the polygons act as subresolution assist features(SRAF) and will themselves not print on a wafer.

Some mask writers are capable of producing patterns having angles otherthan 90°. In this case, an optimized mask data can be approximated usingthe techniques shown in FIG. 24B. Here the additional feature 252 isapproximated by an annular octagon 270 having a number of sidespositioned at 45° to each other. The size of the polygons that make upthe feature 270 act as SRAFs such that they will not print on the wafer.

Line Search Strategy

With various implementations of the invention, the pixel inversion issolved iteratively. More particularly, a first step to identify a searchdirection may be performed, followed by a second “line search” step tothat identifies a point where the objection function becomes close tothe lowest possible value in the local scope of the search direction.

There are various ways to determining the search direction based, suchas for example steepest descent or quasi-Newton. Conventional linesearch strategies may be accomplished by first, scaling the searchdirection (or search vector) by multiplying a predetermined constantscaling factor (called γ) with the search vector. Subsequently, thesearch range may be divided into a preselected number of even sub-steps.Then the objective function may be evaluated starting at the sub-stepclosest to the scaled search vector and ending at the sub-step farthestform the scaled search vector. The first sub-step where the nextsub-step gives a larger objective function value may be selected as theideal candidate.

There are a number of problems with the conventional approaches.Particularly, there is no way of knowing if the search range based on afixed scaling factor will be appropriate for all iterations steps withdifferent search directions, as a result it can lead to unstableconvergence behaviors. Additionally, it is often arbitrary to use afixed integer to divide the search range into the evenly distributedsub-steps. A small sub-step size may find local minima, but it alsotakes more objective function evaluations due to increased sub-steps,thereby making the process slower, whereas a large sub-step may befaster, but it can cause the line search process to miss local minima bya large margin.

FIG. 25 illustrates an objective function ƒ(x+αs) along a line searchdirection. Here, x is a variable vector, s is a search vector, α is aline search scalar parameter, and y is a maximum range of the linesearch function ƒ(x+αs) along a given search direction.

As discussed in [38], the descent property, which is the inner productbetween the search vector and the gradient, may be defined as follows:

sTg(x) where g(x)=δƒ/δx,

which is the first order derivative (or gradient) of the objectivefunction at the starting point of the search vector.

Note: with those algorithms such as steepest descent and quasi-Newton,the first order derivative at the starting point of the search vector iscomputed as part of the overall optimization process. Therefore it canbe assumed that this derivative is available without additionalcalculation during the line search process.

The descent property is the slope of the objective function along thesearch direction at the starting point of the search vector, which isillustrated in FIG. 25B. Based on the descent property, you can get anapproximate objective function value at the search sub-step of as αs asfollows:

ƒ(x+αs)˜ƒ(x)+αsTg(x).

If the expected maximum objective function improvement Δƒ=ƒ(x+γs)−ƒ(x)for the current iteration is known, an estimate as to γ can be made asfollows,

γ=|Δƒ/(sTg(x)ρ)|

where ρ is a slope-relaxation factor, which we will discuss later. Touse this formula to calculate γ, reasonable estimates of Δƒ and ρ shouldbe calculated.

One way to calculate the estimate of Δƒ is to keep track of objectivefunction value improvements in the previous iterations of theoptimization process. The objective function value keeps falling duringthe optimization iterations as illustrated in the “objective functionhistory curve” shown in FIG. 25C.

There is no guarantee that an objective function history curve fallsdown smoothly as the iterations go. It is quite possible the curve fallsrapidly at some points and falls slowly at other points and there is noeasy way to tell what is going to happen in this respect. As a result,estimates of Δƒ could be quite off if the estimate is based on oneparticular iteration result. However, it's generally true that the curveinitially falls rapidly and then later the fall rate significantly slowsdown as it gets closer to convergence. And it's also generally true ifyou take a look at multiple consecutive points as a group, they tend toshow more smoothly falling tendencies. Based on this observation, itshould be a good strategy to make an estimate of Δƒ based on multiplevalues of Δƒ of previous iterations. One simple way to doing this is tocompute the moving average of Δƒ, say for the previous three iterations.

To make an estimate of ρ, one can consider what the value of ρ shouldhave been for the previous iteration. Once a line search process iscompleted for the current iteration of the optimization process, thevalue of real Δƒ can be obtained. The descent property for the startingpoint of the current search vector, sTg(x), is also known.

To see the meaning of ρ in [11], refer back to FIG. 25B, typically theobjective function value along a search direction keeps falling up to apoint, and then after that it starts moving up. As suggested in FIG.25B, the slope of the descent property would be steeper than the slopesat intermediate points in the search direction. To make a betterprediction as to the search range based on some slope like the descentvector, it would be better to make such a slope somewhat less steep asillustrated in FIG. 25D, ρ should be estimated in such a way that themodified slope value of sTg(x)ρ provides a sufficient search range of γfor the expected objective function improvement of Δƒ.

Again, what ρ should have been for the previous iterations can becomputed. To do that, we introduce the value called best α. The best αis the value of α at the point where the objective function value getsthe lowest (or close to the lowest) in a line search process. The best αcan be denoted as A.

If the adjusted slope in FIG. 25D is calculated in such a way that itdirectly points to A, then the slope value σ should be |Δf/A|, asillustrated in FIG. 25E. If you simply use the slope value in FIG. 25Eto compute the search range, the search range is A. But usually you'dlike the search range to be sufficiently larger than the best a so thatpossible points of best α be covered in the search. Say, the searchrange should have been 10 times larger than best α best on the previousiteration result, then

10·A=10·|Δf/σ|.

For the value of (sTg(x)ρ) to cause the 10·A search length based on

γ=|Δf/(sTg(x)ρ)|,(sTg(x)ρ)

should have been equal to |Δf/A|/10. Generally, we put the adjustedvalue of ρ as follows:

ρ=χ|Δf/(sTg(x)A)|

where χ is the aforementioned adjustment factor, which, in the examplementioned above, is equal to 1/10. For similar reasons mentioned beforewith regard to the use of moving average values, we extend the use ofmoving average values to all related parameters. As a result, we canrewrite γ=|Δf/(sTg(x)ρ)| as follows:

γ=|<Δf>/(<sTg(x)><ρ>)|

where <Δf> is the moving average of Δf from the previous iterations, saythe last three, <ρ> is the moving average of p similarly, and <sTg(x)>is the moving average of the descent property from the current iterationplus the previous iterations, say two previous ones.

The actual flow that incorporates this line search control method usingthree point moving averages looks as follows: First, compute the initialobjective function value, ƒ0, and set the initial expected objectivefunction improvement, Δƒ0 equal to ƒ0. Second, start the overalloptimization iterations (e.g. quasi-Newton) with the iteration counterof k (k=0,1,2, . . . ). Third, compute the gradient gk, for the currentiteration. Fourth, compute the search vector sk for the currentiteration. Fifth, compute the current descent property, dk. Sixth, if kis smaller than 3, set γk equal to the following:

|Δƒk−1/dk| if k>0 and k<3

|Δƒ0/dk| if k=0

otherwise set γk equal to the following:

$\gamma_{k} = {{( {( {1/3} ){\sum\limits_{j = 1}^{3}{\Delta \; f_{k - j}}}} )/\lbrack {( {( {1/3} ){\sum\limits_{j = 0}^{2}d_{k - j}}} )( {( {1/3} ){\sum\limits_{j = 1}^{3}\rho_{k - j}}} )} \rbrack}}$if  k ≥ 3

Seventh, run the sub-step evaluation process (described separately) andkeep the following values: Ak, the value of α that gives the lowestobjective function value. ƒk, the objective function value at α=Ak.Eight, compute Δƒk and ρk as follows:

Δƒk =ƒk−1−ƒk

ρk =χ|Δƒk/(sTg(x)Ak)|

Ninth, if the convergence is achieved (e.g. Δƒk has become zero or veryclose to zero) or the maximum number of the iterations is reached, exitthe loop, otherwise increment k and go back to the second step.

As to the sub-step evaluation step (seven), as mentioned earlier, fixedsub-step definition and incremental search from the starting toward theend point may work, but it's not an efficient way if a small sub-stepsize is used, and it's not an accurate way, if a too coarse sub-stepsize is used. One relatively easy improvement to this approach is tostart with a coarse set of sub-steps and identify the sub-step thatcontains the lowest objective function value shown in FIG. 26A. Thendive that sub-step further to conduct the final fine sub-step linesearch FIG. 26B.

Adaptive Weight Adjustment

Once the optimization process is completed, the system has to apply athreshold operation to those mask transmission values to make themattain the discrete mask transmission values that are allowed on a realphysical mask. For example, for the mask whose physical limits are−0.245 (the lowest) and 1.0 (the highest), the system would examine theoptimized mask transmission values on the pixels, and if they are abovethe threshold value of (1−0.245)/2=0.3775 (for example), the finaltransmission values become 1.0 for them, and otherwise they get thevalue of −0.245. The threshold value does not necessarily have to be themiddle point of the upper/lower bounds, and it is indeed possible to usesome other values.

One potential issue in this flow is the possibility that theoptimization process may converge to a state which uses masktransmission values that are not close to the discrete values that areallowed in the final real mask. Such in-between mask transmission valuesare referred to as grayscale values in this document. For example, theoptimization process may decide the mask transmission value of 0.3 isoptimum for some region, and from mathematical point of view, that maywell be a perfectly optimum solution. However, after the thresholdoperation, the pixels in that region will be all brought back to theminimum mask transmission value allowed in the real mask. This could beparticularly problematic for SRAF formation purposes, as there is nosuch thing as grayscale SRAF's in the real mask.

One way to addressing this issue is to add an additional penalty term inthe objective function in such a way those grayscale values cause highercost. However, our experience shows a) whenever you add additional termsthat are quite different in nature from the original optical imagingobjective function, end results tend to become be unreliable andunpredictable from optical optimization point of view, and b) such apenalty term that punishes grayscale values also tend to punish theformation of SRAF's, as those SRAF's tend to be gradually formed takinggrayscale values during the optimization process.

The adaptive weight adjustment method addresses this issue quitedifferently. Instead of adding additional terms to the objectivefunction, this method makes adjustments to the original imaging term ofthe objective function in an adaptive way during the optimization. We'lldiscuss the details in the following section.

As mentioned earlier, the original objective function of the pixelinversion process is defined as follows:

${{f(m)} = {\sum\limits_{ij}{{w_{i,j}( {{I_{i,j}(m)} - T_{i,j}} )}S}}},{{where}\mspace{14mu} ( {{S \geq 2},{{an}{\mspace{11mu} \;}{even}{\mspace{11mu} \;}{integer}}} )}$

When the optimization process converges to a state with lots of grayvalues, that may well be a mathematically valid solution. In otherwords, the objective function as it is simply has not been set up togenerate a desired solution with distinctive SRAF's. Which means, unlessyou modify the objective function in such a way that it leads to adesired solution, there is fundamentally no way that the optimizationprocess could achieve a desired goal. The problem is that it is not atrivial problem to set up the objective function that way. We aredealing with lots of complex geometry shapes in various configurations,not to mention the fact that the whole image is being computed basedadvanced optics/illumination setting. All of those complexities tend toadd up to the situation where it is very difficult to know beforehandhow exactly to set up the objective function to achieve a desiredresult.

There may be a way to solving this problem by adding separate terms tothe objective function. One such additional cost term would be off-tonepenalty, which is defined as follows:

${\sum\limits_{ij}{\lbrack {( {m_{i,j} - m_{\min}} )( {1 - m_{i,j}} )} \rbrack n}},{{where}\mspace{14mu} ( {{n \geq 1},{{an}\mspace{20mu} {integer}}} )}$

Here, m_(min) is the minimum mask transmission value. The problem ofthis approach is a) adding such a term that has nothing to with theoptical behavior could distort the image based pixel inversion results,b) the off-tone penalty tends to keep the transmission values groundedto the two extremity values, thereby making it obstructive to SRAFformation purposes, and c) as a result, end results could beunpredictable and unreliable.

Various implementations of the invention provide a solution to thisproblem which we call the adaptive weight adjustment. This approach isbased on the original formulation of the objective function

${f(m)} = {\sum\limits_{ij}{{w_{i,j}( {{I_{i,j}(m)} - T_{i,j}} )}S}}$

It does not add separate terms to the objective function. However, itdoes change the weight values (wi,j) during the optimization process inan adaptive way.

The basic procedure of the adaptive weight adjustment is as follows:Check the current intensity for each of the pixels, and increase theweight in the following cases:

1. if the intensity for a bright region pixel is lower than the I_(max)by a certain amount, increase the weight for the pixel.2. if the intensity for a dark region pixel is higher than I_(min) by acertain amount, increase the weight for the pixel. Where I_(max) is thetarget intensity for the bright region, and I_(min) the target intensityfor the dark region. As to the selection of the values of I_(max) andI_(min), see the reference [39].

Since we are changing those values of w_(i,j), we are changing theproblem definition itself. However, based on the observation that anill-defined objective function leads to ill-defined results, it makessense to improve the definition of the objective function. And we changethe values of w_(i,j) in such a way that those pixels with undesirableintensity values get penalized more by increasing the values of w_(i,j)for such pixels. This maneuver heightens the objective function value,and creates more room for the optimization process to find a searchdirection that lowers the objective function value for the nextiteration, whereas the optimization process just based on the originalobjective function definition could be stuck in one of the minima.

Adaptive weight adjustment doesn't necessarily directly deal with graybar issues, but a) it should achieve better global convergence withcontinually revised objective function definition, b) it should preventprinting SRAF formation, and c) it is still based on optics only,thereby avoiding unreliability of introducing auxiliary objectivefunction terms.

The best mode formula we use for the weight adjustment is as follows:

If the target intensity for the pixel is I_(max) and the image intensityof the pixel I_(i,j) is smaller than (I_(max)−S), then the new weight iscomputed as follows:

w _(i,j) =w _(i,j) +F·(I _(max) −S−I _(i,j)),

where S (slack) is typically a small positive value, F is a constantfactor with a positive value

Similarly, if the target intensity is I_(min) and the current intensityI_(i,j) is greater than S, then the new weight for the pixel is asfollows:

w _(i,j) =w _(i,j) +F·(I _(i,j) −S).

Otherwise, no weight change is made. There are a couple of parametersused here, namely S and F. Our experience shows the following are areasonable setting for the two parameters:

S=0.01·I _(max)

F may be determined so that the maximum possible improvement in theobjective function due to the adaptive weight change would be about 100in total for the entire pixels.

An example of an optimization flow that incorporates the adaptive weightadjustment is as follows: First, run initial iterations based on theoriginal definition of the objective function. Second, switch to theadaptive weight adjustment mode, and keep adjusting the weights the restof the way. It is important to note, that the timing of the switch tothe adaptive weight adjustment mode is determined by the status aroundthe target edge region. If it reaches the point where there's virtuallyno improvement in the target edge region, then make the switch.Additionally, the objective function value is recomputed with theadjusted weights, and that value is used as if it were the previousiteration's objective function value.

FIGS. 27A and 27B illustrate a 170 nm isolated pattern, with optics:λ=193, NA=0.82, annular, σ_(out)=0.864, σ_(in)=0.54, 6% attenuated PSMbackground. FIG. 27A is without adaptive weight adjustment, while FIG.27B is with adaptive weight adjustment, in which case SRAF's aresuccessfully generated. Whereas in FIG. 27A, the optimization convergesto a state with grayscale values and no real SRAF's are seen after thethreshold operation.

FIGS. 28A and 28B illustrate an 80 nm by 80 nm square contact pattrnwith random placements, where λ=193, NA=1.2 (wet), quasar, σ_(out)=0.93,σ_(in)=0.695, angle =30, 6% attenuated PSM background. FIG. 28A iswithout adaptive weight adjustments while FIG. 28B is with adaptiveweight adjustments. As can be seen from these figures, no SRAF's areprinting when adaptive weight adjustments are used.

Post Pixel Inversion Mask Rule Constraint

Raw optimization results of the pixel inversion results generallycontain polygons with lots of small features and arbitrarily anglededges that break all sorts of mask manufacturing rules. It'sparticularly problematic to have lots of small figures that cause notonly mask-writing inaccuracy but also significant increase in the EBshot count. It is known that about ⅓rd of the mask cost come from the EBwriting time, which is roughly proportional to the EB shot count for themask. Given that, it seems quite prohibitive to have a shot countincrease of 100×, for example.

Geometrical mask simplification has been employed to address this issue[39]. Various implementations of the present invention provide a postpixel inversion MRC process to address this issue.

While the geometrical mask simplification as a post process to the pixelinversion helps improve the manufacturability to some extent, it doesnot necessarily ensure that the result be clean in terms of MRC. Thepost pixel inversion MRC process takes the final pixel inversion resultand tries to turn it into MRC-clean result while ensuring that resultantmask shapes would not cause undesirable imaging effects such as printings-bars.

The purpose of the post pixel inversion MRC process is to transform thefinal result of the pixel inversion into the one with MRC-cleangeometrical shapes while ensuring those modified features do not turnout to cause undesirable effects such as printing SRAF's.

As described in [11], the pixel inversion process involves imageanalysis for the pixels in the work region. We use two schemes for theimage analysis: 1) FFT-based method, and 2) electrical field cachingmethod [11]. The electrical field caching method is suitable for theline-search process of the optimization, and it enables fast evaluationof the objective function. To perform the electrical field cachingmethod, it is necessary to calculate the image of the starting point ofthe line search as well as that of the endpoint of the line search basedon the FFT-based method. Once the results of these two points are known,it is possible to exploit the linearity of the kernel convolution toconduct linear interpolation for the calculation of intermediate stepbetween the two points. Since there is no need for redoing FFT duringthe linear interpolation process, it achieves faster evaluation of theobjective function.

The electrical field caching method can be a useful tool for the purposeof examining the image after the geometrical modification due to MRC.The MRC process itself is a geometrical shape transformation which isnot aware of the impact the transformation causes in terms of imaging.The basic idea of the MRC process is to take the final result of thepixel inversion and the geometrical MRC cleanup results as the startingpoint and the endpoint of the line-search process, respectively. Inother words, compute the difference between the MRC result and the pixelinversion result and regard that as the search vector for the ensuingline search process. Then perform the line search process by using theelectrical field caching method to quickly determine the intermediatestep where the objective function gets larger than that the pixelinversion result. This is a systematic process that determines the pointduring the line search where the MRC cleanup result starts adverselydeviating from the image of the final pixel inversion. We call thisprocess quasi line search.

Based on the idea outlined above, we will describe the baselinealgorithm of this scheme below. In the algorithm description, when theword pixel is used, that means the pixel representation of the mask isused, where mask transmission values could be grayscale and/or discretevalues. When the word geometrical is used, that means the geometricalrepresentation of the mask is sued with physically realizable discretemask transmission values.

In various implementations, the inversion is provided as follows:

1.Start with the final state of the pixel inversion (note: the pixeldata here still has grayscale values).2.Convert the pixel data to the geometrical data based on the specifiedthreshold (typically (1 +m_(min))/2, where m_(min) is the minimum masktransmission value).3.Convert back the geometrical data to the pixels and compute theobjective function value (note: the pixel data here has the discretemask transmission values).4. Run the orthogonal geometrical simplification for the result of 2.5. Perform geometrical MRC cleanup (more about this step is describedlater).6. Convert back the current geometry to pixels.7. Evaluate the objective function for the current state.8. If it's smaller than or nearly the same as the objective functionvalue of the pixel inversion result then exit.9. Otherwise, compute the difference between the MRC's pixel state andthe pixel inversion's pixel state and treat it as a search vector.10. Then run the quasi line-search process using the aforementionedsearch vector and electrical field caching.11. Keep the result right before the objective function value getslarger than the objective function value of the pixel inversion.12. Convert the result to geometrical representation.13. Run the orthogonal geometrical simplification to the result.

As to the geometrical MRC cleanup mentioned in the algorithm, it couldbe implemented in different ways. The baseline approach we employed forthe geometrical MRC cleanup is as follows:

1. Put the polygons that interact with the original target shapes in thecategory of the main feature polygons, and all other polygons in thecategory of SRAF polygons.2. Perform the over-sizing operation by the clearance distance on themain feature polygon to form the clearance area (the clearance distanceis the minimum distance allowed between SRAF polygons and the mainfeature polygons).3. Remove part of the SRAF polygons that are inside the clearance area.4. Combine the remaining SRAF polygons and the main feature polygon asthe current set of polygons.5. Compute the space/width distances for all edges of the current set ofpolygons within the maximum distance of the space/width constrains, andthe results are a set of the pairs of the edges that are within theconstraints distance with one another.6. Examine each of those pairs one by one and compute the distancebetween the two edges.7. If the two violation edge's distance is close enough (by a specifiedamount), then a) if it's for the width violation, then delete the partof the polygon between the two edge from the entire polygon that the twoedge belong to, or b) if it's for the space violation, fill in the gapbetween the two edges to form continuous polygon between the two edges.8. If the two violation edges' distance is not close enough, then a) ifit's for width violation, then move the two edges farther apart tocreate sufficient width between them, or b) if it's for space violation,then move them farther apart to create sufficient space between them.

This is just an example of geometrical MRC algorithm. In reality, it canbe much more complex than this example with more MRC constrains such asarea constraint, etc.

FIG. 29A illustrates a target contact hole (red hatched boxes) andassociated raw pixel inversion result (gray shapes), while FIG. 29Billustrates the target contact hole (red hatched boxes) and associatesorthogonal simplification result (yellow lines) for the raw pixelinversion result. Note the simplification algorithm tries to preservethe raw pixel inversion result including tiny shapes. Additionally, FIG.29C illustrates the target contact hole (red hatched boxes) and the postpixel inversion MRC result (green lines). Note, overall, the MRC resultattains even simpler shapes, and also too tiny shapes are removed.

CONCLUSION

Although certain devices and methods have been described above in termsof the illustrative embodiments, the person of ordinary skill in the artwill recognize that other embodiments, examples, substitutions,modification and alterations are possible. It is intended that thefollowing claims cover such other embodiments, examples, substitutions,modifications and alterations within the spirit and scope of the claims.

APPENDIX

The following references are all incorporated herein by reference, andmay be referred to in the specification above.

-   1. Barouch, E., et al., “Illuminator Optimization for Projection    Printing,” SPIE 3679:697-703, 1999.-   2. Cobb, N., and A. Zakhor, “Fast, Low-Complexity Mask Design,” SPIE    2440:313-327.-   3. Cobb, N.B., “Fast Optical and Process Proximity Correction    Algorithms for Integrated Circuit Manufacturing,” Dissertation,    University of California at Berkeley, 1998.-   4. Erdmann, A., et al., “Towards Automatic Mask and Source    Optimization for Optical Lithography,” SPIE 5377:646-657.-   5. Fienup, J. R., “Phase Retrieval Algorithms: A Comparison,” Appl.    Opt. 21:2758-2769, 1982.-   6. Gamo, H., “Matrix Treatment of Partial Coherence,” Progress in    Optics 3:189-332, Amsterdam, 1963.-   7. Gerchberg, R. W., and W. O. Saxton, “A Practical Algorithm for    the Determination of Phase From Image and Diffraction Plane    Pictures,” Optik 35:237-246, 1972.-   8. Gill, P. E., et al., “Practical Optimization,” Academic Press,    2003.-   9. Golub, G., and C. van Loan, “Matrix Computations,” J. Hopkins    University Press, Baltimore and London, 1996.-   10. Gould, N., “Quadratic Programming. Theory and Methods,” 3^(rd)    FNRC Cycle in Math. Programming, Belgium, 2000.-   11. Granik, Y., “Solving Inverse Problems of Optical    Microlithography,” SPIE, 2005.-   12. Granik, Y., “Source Optimization for Image Fidelity and    Throughput,” JM3:509-522, 2004.-   13. Han, C.-G., et al., “On the Solution of Indefinite Quadratic    Problems Using an Interior-Point Algorithm,” Informatica    3(4):474-496, 1992.-   14. Hansen, P.C., “Rank Deficient and Discrete Ill-Posed Problems,”    SIAM, Philadelphia, 1998.-   15. Hwang, C., et al., “Layer-Specific Illumination for Low k1    Periodic and Semi-Periodic DRAM Cell Patterns: Design Procedure and    Application,” SPIE 5377:947-952.-   16. Inoue, S., et al., “Optimization of Partially Coherent Optical    System for Optical Lithography,” J. Vac. Sci. Technol. B    10(6):3004-3007, 1992.-   17. Jang, S.-H., et al., “Manufacturability Evaluation of    Model-Based OPC Masks,” SPIE 4889:520-529, 2002.-   18. Liu, Y., and A. Zachor, “Binary and Phase-Shifting Image Design    for Optical Lithography,” SPIE 1463:382-399, 1991.-   19. Liu, Y., and A. Zachor, “Optimal Binary Image Design for Optical    Lithography,” SPIE 1264:401-412, 1990.-   20. Luenberger, D., “Linear and Nonlinear Programming,” Kluwer    Academic Publishers, 2003.-   21. Minoux, M., “Mathematical Programming,” Theory and Algorithms,    New York, Wiley, 1986.-   22. Nashold, K., “Image Synthesis—a Means of Producing Superresolved    Binary Images Through Bandlimited Systems,” Dissertation, University    of Wisconsin, Madison, 1987.-   23. Oh, Y.-H., et al., “Optical Proximity Correction of Critical    Layers in DRAW Process of 12 um Minimum Feature Size,” SPIE    4346:1567-1574, 2001.-   24. Oh, Y.-H., et al., “Resolution Enhancement Through Optical    Proximity Correction and Stepper Parameter Optimization for 0.12 um    Mask Pattern,” SPIE 3679:607-613, 1999.-   25. Pati, Y. C., and T. Kailath, “Phase-Shifting Masks for    Microlithography: Automated Design and Mask Requirements,” J. Opt.    Soc. Am. A 11(9):2438-2452, September 1994.-   26. Poonawala, A., and P. Milanfar, “Prewarping Techniques in    Imaging: Applications in Nanotechnology and Biotechnology,” Proc.    SPIE 5674:114-127, 2005.-   27. Qi, L., et. al., “Global Minimization of Normal Quartic    Polynomials Based on Global Descent Directions,” SIAM, J. Optim.    15(1):275-302.-   28. Rosenbluth, A., et al., “Optimum Mask and Source Patterns to    Print a Given Shape,” JM3 1:13-30, 2002.-   29. Saleh, B. E. A., “Optical Bilinear Transformation: General    Properties,” Optica Acta 26(6):777-799, 1979.-   30. Saleh, B. E. A., and K. Nashold, “Image Construction: Optimum    Amplitude and Phase Masks in Photolithography,” Applied Optics    24:1432-1437, 1985.-   31. Sandstrom, T., et. al., “OML: Optical Maskless Lithography for    Economic Design Prototyping and Small-Volume Production,” SPIE    5377:777-797, 2004.-   32. Sayegh, S. I., “Image Restoration and Image Design in Non-Linear    Optical Systems,” Dissertation, University of Wisconsin, Madison,    1982.-   33. Shang, S., et. al., “Simulation-Based SRAF Insertion for 65 nm    Contact Hole Layers,” BACUS, 2005, in print.-   34. Socha, R., et al., “Contact Hole Reticle Optimization by Using    Interference Mapping Lithography (IML),” SPIE 5446:516-534.-   35. Sorensen, D.C., “Newton's Method With a Model Trust Region    Modification,” SIAM, J. Num. Anal. 19:409-426, 1982.-   36. Takajo, H., et. al., “Further Study on the Convergence Property    of the Hybrid Input-Output Algorithm Used for Phase Retrieval,” J.    Opt. Soc. Am, A 16(9):2163-2168, 1999.-   37. Vallishayee, R.R., et al., “Optimization of Stepper Parameters    and Their Influence on OPC,” SPIE 2726:660-665, 1996.-   38. Fletcher, R., “Practical Methods of Optimization,” John Wiley &    Sons, pp. 33-44.-   39. Huang, C. Y., et at., “Model based insertion of assist features    using pixel inversion method: implementation in 65 nm node”, Proc.    SPIE 6283, 62832Y (2006).

1. (canceled)
 2. A computer system including one or more processors thatare configured to execute a sequence of program instructions forcomputing mask data for the creation of one or more photolithographicmasks that print a desired layout pattern on a wafer with aphotolithographic printing system that is designed to print features ata predefined tightest pitch pattern, wherein the instructions cause thecomputer to: read all or a portion of a desired layout pattern; define aset of mask data having a number of pixels that are assigned atransmission value; determine an objective function that compares asimulation of the image intensity on a wafer to an ideal imageintensity; define a set of ideal image intensities from the desiredlayout pattern, wherein the maximum ideal image intensity is selected tobe substantially equal to the maximum image intensity determined from atest image of features at the tightest pitch pattern produced by thephotolithographic imaging system; and minimize the objective function todetermine the transmission values of the pixels in the mask data thatwill produce the desired layout on a wafer.